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THE THEORY OF SYMMETRICAL GRAVITY WAVES 
OF FINITE AMPLITUDE 


STEADY, SYMMETRICAL, PERIODIC WAVES IN 
A CHANNEL OF FINITE DEPTH 


J. GOODY and T. V. DAVIES 
(King’s College, London) 


[Received 30 June 1955.—Revise received 12 April 1956] 


SUMMARY 

The perfect fluid problem of two-dimensional finite amplitude gravity waves in 
a channel of finite depth has been investigated and the four quantities, wave velocity 
(c), wave amplitude (a), wavelength (A), and depth of the liquid (h), are shown to 
depend upon two parameters u and k, where 0 < u < u,,, and 0<k <1. When 
u Umax the wave is on the point of breaking. Tables have been prepared of four 
combinations of c, a, A, and h against uv and k from which it is possible to determine 
(say) the wave velocity when the three quantities a, h, and A are given. An analytical 
formula is also derived for the relation between a, h, and A when the wave is on the 
point of breaking. 


1. Introduction 

Part III (1) of this series of papers on symmetrical periodic gravity waves 
of finite amplitude establishes an approximate solution of the problem for 
a channel of finite depth, but no attempt was made in that paper to 
examine the results numerically. It is proposed in this paper to give tables 
and diagrams from which the wave velocity can be derived when the 


height and the length of the wave are known together with the depth of the 


channel and the liquid speed at one point. 

The theoretical problem has received attention in the past from Korteweg 
and De Vries (2) and from Lamb (3), who uses the Rayleigh—Boussinesq 
method, and has developed what has been called the cnoidal wave (since it 
involves the Jacobi elliptic function cn); the solution of Part IIT is also 
expressed in terms of Jacobian elliptic functions, but not only is the nature 
of the approximation clearer in this investigation but it provides additional 
information, in that it is possible to discuss with relative simplicity the 
limiting or breaking wave. As far as the authors are aware no numerical 
results have been given for the cnoidal waves and no detailed comparison 
is therefore possible. It is clear, however, that the present solution bears 

+ Parts I, II, and III have been published elsewhere. 
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the same relation to the cnoidal solution as Packham’s solution (4) of the 
solitary wave problem bears to the Rayleigh solution (cf. 3). Comparison 
with earlier work is possible in the two special cases when the channel is 
infinite in depth and the wave is infinite in length. Parts I and II have 
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dealt adequately with these comparisons; the only results which need be 
mentioned are that the first approximate solution obtained by the present 
technique over-estimates the wave velocity and gives also an over-estimate 
of the critical ratio a/A for the infinite depth problem and of the critical 
ratio a/h for the breaking solitary wave (see Fig. 1). 

The axis Ox is horizontal, Oy is vertically downwards, and the origin O 
is always at a crest moving with a velocity c from right to left (see Fig. 1). 
With respect to these axes the motion is steady. 

A velocity potential ¢ and stream function % exist, and are defined by 


dd = u* dx+v* dy, 
dis -v* dx+tu* dy, (1.1) 


where (u*, v*) are the velocity components relative to O(2, y). 

By symmetry, OF, BC, AD,,... are lines of constant ¢; let 6 = 0 on OF, 
¢ = d, on BC, so that 6 = —d, on AD. Also let the bed CFD and the 
free surface AO B, which are streamlines, bef = ys, andy = 0 respectively. 
The wavelength AB is A, the depth of the bed below a trough is BC = h, 
the amplitude OF = a, and the velocity of the liquid at F, relative to a 
fixed frame, is U. 
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THEORY OF SYMMETRICAL GRAVITY WAVES 3 
e | The problem to be investigated can be stated as follows: knowing the 
n values of U, a, h, and A, to find c. 


It is possible to establish two formulae of the form 
E,(a, h) E(u,k), F(a, A) F(u, k), (1.2) 
together with one of the form 
D,(a, U,c) D(u, k), (1.3) 


where wv and & are two parameters; thus when U, a, h, and A are given, 


w and & can be determined from (1.2) and we may calculate c by (1.3). 


2. The complex potential 


From (1.1) we have 


dw (u* —iv*) dz. (2.1) 
If q (= ce’) is the velocity of a fluid particle relative to O(x, y), and @ is the 


inclination to the a-axis, then 


uv 


7+10. (2.2) 


Let the velocity at O be gy, then Bernoulli’s equation on the free surface 


gives 


q° de | 2qy, us 0: (2.3) 
be this may be transformed, following Levi-Civita (5), into the form 
nt 
co | mg 
ite —e-* sin G, us 0. 
Cr 3 
al 
In Part III this is replaced by the approximate condition 
() 
oo gl ‘ — 
I ¢-37 sin 36, = 8, (2.4) 
] Oub 3 
; where / takes a suitable numerical value, or by 
»\ 
dé gl ae " 
im) ——_ J est 0, ub 0. (2.5) 
| dw c3 | 
At the crest and on the bed of the channel the velocity is in a horizontal 
direction ; thus we are required in this problem to find an analytic function é 
of win 0 < & < &, having a real period 24, and satisfying the following 
= 
I conditions 
he 
(1 imé 0 when @ 0 or + db : 
ly 0 


(11) imé 0 when ws wy: 
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It is shown in Part III that the solution is given by 


f= az) ( , ait —k* on( 4, Asn FG dabs,), k r (2.6) 


where sn(u,k) is one of the Jacobian elliptic functions with modulus 4, 
4K is the real period of these functions and, in order to satisfy (2.4), it is 
necessary that the following relation should hold between the parameters 


3lgys, 2u dn(2u, k’) 


(U+c) sn(2u, k’)en(2u, k’) 


A(u, k), (2.7) 


where u = Ky,/¢9. In (2.7) the modulus k’ is related to k by 

e442? = | (2.8) 
and 21K’ is the so-called imaginary period for the Jacobian elliptic function 
sn(u,k). Furthermore, the wave exists in the range 0 < wu < 4K’. When 
u is small the wave amplitude to wavelength ratio is small, and when 
u 1K’ the wave is on the point of breaking at the crest, possessing at 
this point an angle of 120 degrees. 


3. Formulae obtained by integration of dw/dz 
We may rewrite (2.6) in the form 

Iz } 

( ; | 


Ute = 
( a | 


»[2i Ky, \-2 . 
k2sn2{= t)on? (w— trb,) (3.1) 
rz do 


each elliptic function having modulus k. If we integrate this along OF in 
Fig. 1 we obtain 
oh 


(U-+c)(a+h) = a {1 —22sn2 2( 2a, kone 


0 


io—¢,), k| i diy. 


0 


where, as earlier, u = Ky,/¢). By a change of variable, and by applying 
Jacobi’s imaginary transformation, this relation can be written in the form 


( U+ c)(a + h ) 


10 
= =) | i -sn2(2u, k’ eae os Ma € = B(u,k). (3.2) 
0 
Another formula can be obtained by integrating dz/dw along the bed FC 
in Fig. 1 and we obtain 
$0 
k\(U +c) = | {1—k*® sn2(2iu, k)sn®( K/h, k)}-* dd. 


0 
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3y a change of variable and by applying again Jacobi’s imaginary trans- 
formation this relation can be written in the form 
K 
; (U+c)A eni(2u,k’) ‘ fate Loe ' 
C, nr {1—sn?(2u, k’)dn?(&, k)}-? dé = C(u,k). 
mY) ( 
0 (3.3) 
A fourth independent result follows from an integration along the free 
surface in which an integral expression may be obtained for the amplitude a 
ifthe wave. It is more convenient, however, to derive this result approxi- 
mately, as in (1), by using Bernoulli’s equation on the free surface (2.3). 
From (2.6) we obtain along the free surface 


dw : . 
; (U+e)1—k sn*(27u, k)sn” 
a2 | 





K . 4 
(¢d i) | 


and thus at O, the crest, and B, the trough, we have 
Yq = (U+e){1—k* sn?(2iu, k)sn(iu, k)}, 
dp = (U+e){1—k* sn?(2tu, k)sn?(K —iu, hk)! 
By applying the addition formula and the Jacobi imaginary transformation 
the above velocities may be written 
2 sn*(2u, k (uw, k’)\3 
do (l c)\1 Pri 
en?(2u, P) (u, k’)|’ 


k*sn?(2u,k’) — \} 
1B (U+ex t+ —. eee na 
| en?(2u, k’)dn?(u, k’)} 


If we specialize (2.3) to the point B we obtain 


D 2ga I k? sn?(2u, k’) \5 
ss (U+e) | en?(2u, k’)dn?(u, k’)| 
()_ pe sn?(2u, k’)sn?(u, k’))§ D(u, k) (3.4) 
~ en?(2u, k’)en?(u, k’)| i 


The four functions A, B, C, and D in (2.7), (3.2), (3.3), and (3.4) have been 


calculated for values of 7 and 0, where 


‘Kk’ 
u é ; k cos 6 
L8O 
7) r(10) < 60. 0° < O15") < WO’, 


and the results are tabulated in the next section. It may be mentioned that 
the case r = 60, or u = $k’, is that of the limiting wave which is on the 
point of breaking at the crest. The limiting or breaking wave is of sufficient 
interest to warrant an analytical investigation and the present section will 
be concluded with a description of the breaking wave whose wavelength is 
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long compared with the depth of the liquid. In the case of infinite wave- 


length the value of the ratio h:a for the limiting wave has been given by 


Packham to be 1-21 (4) and by McCowan to be 1-28 (6), and it is proposed 
here to derive the generalization of this property for large 4, which corre- 
sponds here to values of k which are near unity or to values of k’ which 
are near zero. Accordingly it is necessary to obtain the expressions for A, 
B, C, D for small values of k’. 








When u = 4K’ and k’ is small, the following results are obtained from 
2.7), (3.2), (3.3), and (3.4) 
2 9k’ 2k’ J! as 
Sg re a me ION; (8.5) 
(U+c)® 3sn(3K’,k’)en(3K’',k’)  3v3 
K’ 
ee. { 8 eni(2K’ k’) cr : 2¢ fe’)) 1 
(U+e)(a A). mee ©) | I sn?(3K’, k’) ane &)) 3 dé 
wb, K J | en?(é, k’)| 
0 
3.2! 
3.2 hy + O(k's) (3.6) 
$77 
where a ( (1—?sec?a)-} da; (3.7) 
0 


kK 
(l | c)A 3oni(gk", k’) | f] sn2(2K’, k’)dn2(€, k) i dé 





2, kK’ J 
0 
3.23 k’ 
- ~—. jog | —| -+- 1-7306: 3.8 
= oe (3) (5.8) 
2ga { ios k? sn?(2K’, k’) \5 
(U+c)? | en®(2K’, k’)dn2(2K’, k’)| 
4i{] +144 O(k'4)). (3.9) 


In deriving (3.8) considerable use has been made of some integrals given 
by Greenhill (7). Choosing / = } in (3.5), as in Part I (8) of this series, 
it follows from (3.5), (3.6), and (3.9) that 


a+h 2k, 


VLLIk2 LOD (3.10) 
2a V3{1-+ $k’ + O(k'4)} 
Similarly from (3.6) and (3.8) it follows that 
2(a+} 
e v) 2 f -log 4h’ | 1-438!-1, (3.11) 


and thus by eliminating k’ it follows that the lengths a, h, and A are related, 
in a breaking wave, by the formula 
h 4], 


. oe eo (3.12) 
a V3{1+35-48 exp[—J,A/(a+h)]} 





Usin; 
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Using numerical integration we find that J, = 1-066 and thus 
I 2-46 
l4— = : —n (3.13) 


a 1+35-48exp[—1-066A/(a+h)] 
\s A > 0 this gives h/a 1-46 which must be compared with the values 
-21 and 1-28 due to Packham (4) and McCowan (6) respectively. The 
difference in the value of this ratio will be discussed shortly, but it is 
important here to emphasize that the effect of a finite wavelength is to 
decrease the critical value of the ratio h/a. Since h/a > 0 it is necessary 
that A/(a+h) > 2-99 for the above formula to apply. 

In Packham’s derivation of the critical ratio h/a = 1-21, the method of 
procedure is slightly different from that given above. Formulae (3.5) and 
3.9) are used by Packham but, in place of the formula (3.6) a more exact 
formula is used which is derived from the condition that there is no absolute 
motion at infinity; this implies that 4, — ch and U satisfies 


(U+c)® = c* cos? jx. 
Using these two results together with 
gis, 4a 2ga 
(U+c)?  3v3’ (U+c)? 
t follows that ha |-21. Similarly, in the present problem, we could 


postulate that the absolute velocity at a point on the horizontal bed 

beneath a trough is zero, but there would be no justification for such a 

procedure 

4. The determination of the velocity of propagation of the waves, 
given the physical dimensions 


We have already introduced the four functions 





A ek et 
(( | cs 
B, << Shes, = B(u,k) 
wy 
wie (4.1) 
C., ( +-¢) C(u, k) 
Ds, 
a _ Diu, k 
D), ahi ee (u, k) 
and we define two new functions E and F as follows: 
; AB 3l(a+h) : 
T(u,k) = . E 
Kk ) D . oJ 
10 3A . ats 
F(u, k) -- : - == : F. 
D 4da ’ 
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Knowing uv and k, A and D are easily evaluated numerically using tables 
of elliptic functions, etc. The integral that occurs in B can be found by 
Simpson’s rule in all cases except 7 = 60 or u = 4K’, which is the breaking 
case. In this latter case the integrand is infinite at € = wu and the integral, 
which is convergent, has been found by expanding the integrand in ascend- 
ing powers of (—w) over some suitable range u—e to wu and using Simpson's 
rule over the remainder of the range. It is most suitable in the integral of 
C to expand the integrand. For k > } it can be expanded as a power series 
in dn*(é, /), each term being integrable by a reduction formula. For k < } 
the integrand can be expanded as a fairly rapidly converging power series 
in k*. 

Having found A, B, C’, D for appropriate ranges of values of wu and k, 
E and F can be found. The values of # and F for given values of u and k 
have been tabulated and graphs of H = constant and F = constant have 
been drawn in the (7, @)-plane. 

In a particular problem, if we know the physical quantities a, h, A, the 
values of # and F, say E, and F, respectively, can be found from (4.2). The 
graphs of EH = E, and F = F, are then superposed on the same (7, @)-plane 
and their point of intersection, if any, determines the values of the para- 
meters r and @, or alternatively u and k for the problem. The value of D, 
that is (U’+c), may then be read from the tables. 

TABLE 1 
Values of the function A 


























6 
y ° 15 30 $5 60 75 go 
1'0000 1*0000 I*0000 1*0000 1*0000 1*0000 
10 1'0206 1°0472 1'0207 1°0215 1'O241 1°0337 
20 1:0861 1'0861 1°0867 170895 1°0994 1°1352 
30 1*2092 1*2092 1*2104 1°2160 1°2361 1°3083 
40 1°4178 1°4179 1°4196 1°4281 1°4586 1°5679 
50 1°7723 I°772 1°7745 1°7551 1°8232 1°9599 
0 2°4184 2°4185 2°4207 2°4320 3°4722 2°6166 
a al — © 
TABLE 2 
Values of the function B 
6 


60 











e 


THEORY OF SYMMETRICAL GRAVITY WAVES 


TABLE 3 


Values of the function C 





30 45 690 75 ge 
22°9397 17°9303 
1°3331 88559 
5 7° 3998 5°7725 
} 5°3452 $°1740 } 
+ fOIS7T 3°1527 
106 2° 3966 








TABLE 4 


Values of the function D 








45 60 75 ge 
“0000 "0000 0*0000 0*0000 0*0000 
0181 0°01 47 00100 0°0046 0*0000 
8 0782 0:0648 0°0528 0°0229 0*0000 
75 2036 O1731 O*1 301 0°07 34 0°0000 
4 } *4403 0°3951 0°3175 0°2082 00000 
- <Sy 0°8772 7574 0°5807 0*0000 
8 { 4332 2°3238 2°1674 1°9543 00000 
TABLE 5 
Values of E and log,, E 
A 
45 60 75 )O 
7 | 4035 69°4997 224°7174 
7413 1°8420 2°3516 
"9034 16°8217 20°8303 49°5721 
1421 1°2261 1°3187 1°6953 
{ 657 7°O424 5239 17°8563 
: 7757 0°3477 0°9788 1°2518 
47 1999 yO514 4°6350 7°5925 
4 4754 sOsl 0°5625 06604 08804 | 
8 )210 2°1105 2°4294 3°4537 
2825 0°3244 O°2355 0°5420 
I 1°3165 1°4320 1°55 
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TABLE 6 


Values of F and logy, F 


> 








0 
° [5 30 45 60 75 90 

fe 1625°34 1293°62 1245°97 1437°65 2332°72 
3°2109 31119 370955 3°1576 3°3679 

20 are I192°20 157°49 148-90 144°96 256°33 
2°2535 2°1973 2°1730 2°1612 2°4055 

30 52°5377 4$3°9755 40°5509 42°7252 61°0213 
1°7205 1°6432 1°6080 1°6307 19855 

jo 20°5312 16°9263 15*0570 15°3379 19°0347 
1°3124 1°2285 | 11785 1°1857 1°2795 

50 5°5590 7°4430 | 6°4157 61941 6°6313 
0°9324 o'8717 0°8073 0°7920 0°8216 

5o 3°7202 | 370535 2°5072 2°2537 2°0765 
0°5706 074848 0°3992 0° 3529 0°3173 








Applications and examples illustrating the use of the tables 

In a particular case it may be found that the appropriate HZ and F curves 
do not meet : in this case the progressive wave of that description does not 
exist. Take, for example, 


h A 
= ()-2, 12-6; 
a a 
we then have #, = 1-2, F, = 3-16 (in numerical examples the value of / is 
taken to be 4) and thus 
Ly” 1 v 1 
logy Ly = is, logyo fy = 3- 


[t will be found that these two curves do not meet and so there is no pro- 
gressive wave with these dimensions. 

It should be noted that the range of values of E and F give minimum 
values of h/a and A/a. The minimum value of E is approximately 1-2 and 
so h > !a; hence using (3.13) 


0-68h < a < 5h. 


Similarly, from the tables, the minimum value of F is greater than 2, hence 


A > 8a. This result must be compared with the corresponding one in 
Part I (8), namely a/A = 0-126. 
Consider now the case 


Aja = 40, hia 





the 


Th 
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2:0 
10 
4041| \ 
0:8 
. 
" 06 
0-4 
7 0:3 
‘ 302C «O45 (itiwsti«<CTSts«SD 
A 
Fic. 2. Graphs of curves 1/log,, F J for J 0-3, 0-4, 0-6, 0-8, 1-0, 2-0. 


10 











1! cn: CL nn 10) 
—A 
Fic. 3. Gs iphs of curves | OL19 BE J for J 0-5, O-7, 1-0, 1-5, 2-0, 4-0, 8-0. 

then £, 1-6, Fi 10 and thus 

| 7 l 

: 1-5, - he 
logs E, logio i 

These curves meet at the point given by 0 = 75°, r = 48. Then, either 


by interpolation between the tabulated values, or by direct calculation, 
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we find that D = 0-525. Hence the velocity of propagation is given by 


the formula > 
2ga 


(U+c)? © 


0-525. 


The velocity of propagation may be determined only when a and U’ are 
given. If the wave amplitude is 1 inch, so that a =}, , we have 

U'+e = 3-2 feet per second. 
The fact that c is indeterminate as long as U is unknown in this finite depth 
problem was pointed out originally by Stokes, and it is clear that this 
difficulty will remain as long as the liquid is assumed to be non-viscous. 
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FLOW WITH VARIABLE SHEAR PAST 


CIRCULAR CYLINDERS 
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SUMMARY 


Stream functions are obtained for two variable shear flows past a circular cylinder. 
The shear distributions in the free stream are of the ‘linear’ and ‘boundary layer’ 
type respectively. In both flows, the stagnation streamline is displaced in the free 
stream towards a region of higher velocity. This is in agreement with the results 
obtained by the present authors (1) in the case of constant shear flow past a cylinder 
and with the experimental results of Young and Maas (2) for a pitot tube of circular 
cross-section in a low-speed shear flow. The displacements in both variable shear 
flows are comparable in magnitude to the values obtained in the constant shear case 


but are considerably smaller than the experimental values of Young and Maas. 


List of notation 


J Bessel function of the first kind 

| modified Bessel function of the first kind 

if Bessel function of the second kind 

K modified Bessel function of the second kind 
l typical linear dimension 

L perimeter of cylinder 

N rotational flow parameter (= U/w,l) 

rey cartesian coordinates 


7,4, Velocity components in the x- and y-directions 
r.@ polar coordinates 


qe radial and transverse velocity components 


q velocity 
If standard velocity 
y. deflexion of the stagnation streamline in the free stream 


circulation 
7) thickness of the incompressible boundary layer 
stream function 
7 disturbance stream function 
Ww vorticity 


Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 1 (1957)] 











14 J. D. MURRAY AND A. R. MITCHELL 


1. Introduction 


THE stream function / for the steady rotational two-dimensional flow 


of an incompressible inviscid fluid satisfies the equation 


Oy Orb 


= ——— —— 0, (1) 
ox" oy” 
where the vorticity w, given by 
cq, 
ne (2) 
ox = ey 
2 : Cw CW > 
satisfies the equation q, qy- 0. (3) 
Ga "oy 
ous ous 
and = —, q, . oY 
oy ox 


are the velocity components parallel to the x- and y-axes respectively. 

As far as the present authors are aware, no exact solutions of (1) and 
(3) exist except in the case of constant rotation, although several investi- 
gations using approximate methods have been carried out. In the present 
paper, the stream functions for some variable shear flows past a circular 
cylinder are obtained, and certain aspects of these flows compared with 
experimental results. 


2. The ‘linear’ shear distribution 


The general solution of (3) is 


w Ff (), 
and so the vorticity is constant along a streamline. In the particular case 
, ub 
F(t) 2? 
"0 


where 7, is an arbitrary length, (1) becomes 


Ons r Cus ub 


AS ee =e (4) 
dat  oy® =r 
A vorticity distribution in the free stream near x —oo, consistent with 
(4), is 
Ww eulro e~UiTo (5) 
9 
rO 9 
which on integration gives the velocity distribution 
a , 
Vr eiro_ e-u ro. (6) 
ro 9 


where the constants a and / are arbitrary. If a standard rotation w, and 
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a standard velocity U are introduced as the values of w and q, where y = 0, 
it follows from (5) and (6) that 


1 


a 37o(U —arg1o), b b7o(U + WoT). 


If these values are substituted in (5) and (6), the vorticity and velocity 


distributions in the free stream near x 1 become 
w y — 1 - 
cosh ~ — N sinh” ; (7) 
Wy ro ro 
q y ere 
and 7 cosh ~— — sinh~, (8) 
ry A "0 
respectively, where .V (’/(we%o). In order to keep the velocity positive 


everywhere the parameter NV must satisfy the inequality N > 1. The 
above distributions are now expanded in powers of y/r,. If higher powers 


are neglected, the expansions become 
Ww x(¥), (9) 
Wy r9 
, l l 2 
and fa ] (?) 4 3(¢] . (10) 
( N r9 =\"o 


respectively, and so the free stream shear distribution (7) is approximately 


ot y/? 


linear in the vicinity of the x-axis. In Fig. 1 the similarity of the velocity 
profiles given by (8) and (10) in the vicinity of the x-axis is shown for the 
case NV 2. 

A circular cylinder of radius 7, is now introduced at the origin and it is 
required to find a stream function which satisfies (4) and takes a constant 
value on the cylinder. For convenience, introduce the function ‘’ where 


b = aevitot be-virot (11) 


and from (4), ’, which tends to zero as x tends to —, satisfies the equation 
er ey y 


» 9 
Cx cy~ 


; (12) 
ro 

Using polar coordinates r,@, and the method of separation of variables, 
12) is solved to give 


7 B, OK (r')+ > (A,,cosmé-+ B,, sin m6)K,,(r’), (13) 


0 
where A,,, B, are constants, A,, is the modified Bessel function of the 
second kind of order m, and r’ = r/ry. Thus from (11) the required stream 


function us is given by 
us aero be YiT¢ B, 6K,(r’) T S (A 
os 


0 


_cos mé+- B,, sin mé)K,,(r’), 


” 
m 
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where y’ = y/ry. Now & must take a constant value, say c, on the circular 
cylinder r = 7), and so 
x 
:; sin® sin 8 , . _ r , 
ce = aesin? + be-sin’__ Bo AK (1)+ > (A,,cosmé+ B,, sinmé)K,,(1). (14) 
m=0 
The Fourier expansion of +s"? is 
x. 


ow 
sin 0 ‘ a ¢ , Je ng 
exsin? — J(1)+2 > (—1)"d5,,(1)cos 2mO+2 ¥ (—1)",,,,(1)sin(2m+ 1)6. 
m=1 m=0 
The above values are substituted into (14) and the coefficients of corre- 
sponding sine and cosine terms equated. With the values of a and } found 
previously, the stream function becomes 


us . ’ l ; Lae ; ] K, (r’) 
¥ — sinhy’——coshy’+|—"—4 1)| 20") 
U;. inh y v° Ly Ee 5 Ih ( } K(1) * 
ahs 2 S (—1)™ Lam 2 K,,,(r’ eos 2mé 
WV ~ K,,,(1) 2m alias 
a. er | anes ‘ 
+2 > (—1) Katy ame” )sin(2m+1)6@. (15) 
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The velocity components in polar coordinates are 
l cds Cds 


oe and qd ——., 
r 00 te or 


Round the circumference L of the circular cylinder, g, = 0, and the circu- 
lation [in an anticlockwise direction is given by 
I D dor, dO +O. 
L 


where [ 


0 


and IT, are the circulations arising from the undisturbed flow and 


the disturbance flow respectively. Using (15), 


l) l 
1,1), 
27rUr, N i) 
> Ko) fp 1 
and * of -[(1) T ? 
2rUr, Ky(1) | A Urs 
where K,(1) (dK, (r')/dr'),._,.. Now in the case of constant shear flow 


past a circular cylinder, Tsien (3) (see also Richardson (4)) uses the hypothe- 
sis that the introduction of a cylinder into an inviscid incompressible rota- 
tional flow does not alter the circulation round the path L, which coincides 
with the perimeter of the cylinder. In the present problem this leads to 


the result [° 0 and so 


l 


( Wo rz I,( 1). 


If this value is substituted into (15), the stream function becomes 





al 
a , ae £. . ’ 
J sinh 7’ - cosh y 6 \ (—1)" am(1) ky, (7 jcos 2mé- 
Ur, 7 HN Wa in 
14 m=1 om 
‘ i i one : 
2> (-1p*=S> al’) Ky ,.4(7")sin(2m+1)@, (16) 
mga Konii1(1) i 
)0 resulting in a velocity on the cylinder of 
: ; = , 
re (Va)w=4 sin 6! cosh(sin @) - sinh(sin 9)| | 
| N 
nt 
: YS (-1)" Fam| I) K4,,(1)cos 2mé - 
N -- Ky,,(1) 
. i i ‘ = 
ya ] m+ Lom +a K... {()sin(2m-+1)8], (17) 
oper Kom ii) i 
where 
dk,,,() > LK, i 
Kot?) Bem! | Kom 4i(l) = am +i(T ) 
. d) | 7 dr 
15) ? 1 " 1 


and gg = q/U. The stagnation points on the cylinder are given by 


092 .37 Cc 
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(9@)-1 = 0, and the stagnation streamline, given by %/Ur, = —J,(1)/N, 
starts in the free stream near x = —o at the position y, given by 
‘ 72(1)+(N?2—1)}f}—L() . fs 
yf. log, | o( ) ( 7 )] of ) (N ss 1). (18) 
ie N—1 
= a. 3 T T — 7 
| | SpHene wm 
| Cowstant SHEAR 
ob } = = 
Youne AND Mans. 
ly’ Shear 
IY. a 
Linear SHEAR 
2 ae t 4 ae 
| CowstanT SHEAR 
Boun oaRY | 
iLayver Swear | 
: | | | 
0 1-0 
0 25 WY 0-5 0-75 
Fia. 2 


The deflexion of the stagnation streamline is shown in Fig. 2. It should 
be noted that the convergence of the series in (16) and (17) is extremely 
rapid and the series may be terminated at m= 4 for all practical 
purposes. 


3. The ‘boundary layer’ shear distribution 
Another particular solution of (3) is 


us 
WwW > 
2 
where / is an arbitrary length. On substitution into (1) this leads to the 
equation ng ai 
juat Ow emus yb (19) 
Ox? by? [2° 
Vorticity and velocity distributions in the free stream near x = —0o con- 
sistent with (19) are , 
a y\ b. fy 
w = —cos|{~]+—sin[=}, (20) 
? ct} P l 


» 


1 
and q,. = ; cos (7) sin (7); (: 
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respectively, where a and } are arbitrary constants. If a standard rotation 
w, and a standard velocity U are introduced as the values of —w and q,, 
where y = 0 and y = 5—y, respectively, it follows from (20) and (21) that 














3) 
— 
P - } N’—sin{(d—y»)/U} 2 
Wl"; ) - a a Wo; 
cos{(d—y,)/L} 
where V’ = U/(wol), and yy) is a variable parameter lying in the range 
) J 5. It 
" Yo Yo ] 
COS —~». 
l A 
_ 
we — 
a 
j 
il i ai 
Fic. 3 
and the arbitrary length / is chosen to be 
ae 
l oO, 
e the vorticity and velocity distributions in the free stream near x = —oo 
become (2 ’ )\ 
) . N’cos(” YT Yo)\ (22) 
Wy 20 
7(Y + Yo) 
and : sin| ——_—}, (23) 
20 
)) 
respectively. In the range 0 < y+y, < 4, the free stream velocity distri- 
bution (23) is Lamb’s approximation for the flow in an incompressible 


boundary layer of thickness 6 (5). 
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A circular cylinder of radius r, is now introduced at the origin, and the 
stream function required which satisfies (19) and takes a constant value 
on the cylinder. Introduce ‘¥ where 


¢ = acos* +b sin + . (24) 


then from (19), ‘’, which tends to zero as 2 tends to —, satisfies the 
equation cour -n0ur : 
| er ey y 


ene 25) 
Ox? Oy? - ( 


An appropriate solution of (25) is 
oe 
q By 0¥)(r')+ > (A,, cosm0+ B,, sin mO)¥,,(r’), 
m=0 
where A,,, B,, are constants, Y is the Bessel function of the second kind, 
and r’ = r/l. Thus from (24) the required stream function ¢ is given by 
Y = acosy’+bsiny’+ B, 6Y,(r’)+ > (A,, cosm6+ B,, sin m8)Y,,(r’), 
m=0 
and since ¢% takes a constant value, say c, on the circular cylinder r = 15, it 
follows that 
ce = acos(ksin 6)+-bsin(k sin 6)+ B, 6Y,(1)-+ 


+ > (A,, cos m6+ B,, sin mé@)Y,,(1), (26) 
m=0 
where k = zr,/26. The Fourier expansions of cos(k sin #) and sin(k sin @) are 
cos(ksin@) = A(k)+2 > Jy,,(k)cos 2mé 
1 


m 


14y(k)sin(2m-+ 1)8, 


" 


f 
and sin(ksin@) = 2 > J, 
m=0 
respectively, where J is the Bessel function of the first kind. The above 
series are substituted into (26) and the coefficients of corresponding sine 
and cosine terms equated. Following the method of the previous section, 
and with the present values of a and b, the stream function becomes 


ip < Jom (k) 


a cos(y’ +45) +2 cos y4 Y,,, (1 cos 2mé 


-2sin yj Ss Fe Yan sa sin +1)0, (27) 
+1 


where y’ = y/l, yo = yo/l, and the constant value taken by ¢ on the cylinder 


1S 


c= —w,l2J,(k). (28) 
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The stagnation points on the cylinder are again given by (é@’/ér’),,_;, = 0 








with w’ u/Ul, and the stagnation streamline starts in the free stream 
near x © at the position Y. given by 
, J, (k l 
y cos~} o(*) cos~! —1, (29) 
N N 
where 7, y./l 


The field of flow represented by the stream function (27) is a good 
\pproximation to the flow past a circular cylinder situated in a two- 
dimensional boundary layer provided r,/6 is small and y,/5 has a value not 
too near the end values of its allowable range 0 < y,/5 < 1. The problem 
is illustrated in Fig. 3. The deflexion of the stagnation streamline is shown 
in Fig. 2 for 7/6 5 (k = 7/40) over the range | < y,/6 < 4, where 


N | U) Wo 1) is equal to (40/77) N’. 


4. Comparison with experimental results 

In measuring the profile drag of a wing-section by the pitot traverse 
method, there is a considerable velocity gradient across the wake traversed 
by the pitot and so it cannot be assumed that the pressure registered by 
the tube is the total pressure of the stream opposite the geometric centre 
of the face of the pitot. Accordingly, Young and Maas (2) carried out tests 
which showed that when a pitot tube of circular cross-section is placed at 
zero incidence in a low speed flow with a constant total pressure gradient, 
the stagnation streamline comes from a region where the velocity is higher 
than the velocity upstream on the pitot axis. Ifd is the upstream displace- 
ment of the stagnation streamline and 7, is the outside radius of the pitot 
tube, Young and Maas found that d/r, is approximately constant for a 
variety of conditions and equal to 0-36. In terms of N defined by 


N = Ulery 


the range covered by the experiments was 0 < 1/N < 0-3, where U and 
vy are the values of the velocity and rotation at the geometric centre of 
the face of the pitot tube. 

This displacement of the stagnation streamline when a pitot is placed in a 
flow across which there is a constant total pressure gradient is shown in Fig. 2, 
together with the magnitudes of the displacements previously obtained 
from theoretical considerations of inviscid incompressible shear flows past 
a circular cylinder. The types of free stream shear distribution considered 
ire constant, ‘linear’, and ‘boundary layer’ respectively. In the constant 


shear case (1), with the usual notation, the displacement of the stagnation 


streamline is given by 


N—(N?+1)}, (30) 
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where V > 0, and U is the free stream velocity at y = 0. Also shown in 
Fig. 2 is the displacement obtained using theoretical methods by Hall (6) 
for constant shear flow past a sphere. 

In drawing conclusions from Fig. 2 it must be kept in mind that in Young 
and Maas’s experiments the shear distribution in the wake, across which 
the total pressure gradient is constant, does not correspond directly to any 
of the theoretical models, although it is closest to the ‘boundary layer’ 
distribution of shear. Also, in the theoretical models the fluid is assumed 
to be inviscid, incompressible, and two-dimensional. Despite the above 
differences between theory and experiment, the following broad conclusions 
can be reached by comparing the experimental with the theoretical results 
in Fig. 2: 


(1) The contention of Young and Maas that the stagnation streamline 
is displaced towards a region of higher velocity is supported by all the 
theoretical results. 

(2) It seems probable that the curve obtained from experiment by Young 
and Maas should pass through the origin of the coordinate system. 

(3) In view of the small range of 1/.V (0 < 1/N < 0-3) covered by Young 
and Maas’s experiments it seems rather premature to assume that d/r, will 
equal 0-36 at values of 1/N considerably in excess of 0-3. Recently, however, 
MacMillan (7) using a circular pitot has traversed the turbulent boundary 
layers adjacent to a circular pipe and a flat plate. For a variety of positions 
of the pitot, he found that d/r, was substantially constant and never in 
excess of 0-32. 

(4) The stagnation streamline displacements obtained from the ‘linear’ 
and ‘boundary layer’ type shear flows past a circular cylinder are com- 
parable in magnitude to the displacement obtained in the constant shear 
case. This suggests that in two-dimensional shear flow problems the magni- 
tude of the displacement is not greatly influenced by the type of shear 
distribution in the free stream. 

(5) At small values of 1/N the experimental displacement of Young and 
Maas is much greater than the displacements obtained from the theoretical 
solutions of two-dimensional shear flow problems. It is, however, com- 
parable in magnitude to the theoretical value obtained by Hall for constant 
shear flow past a sphere, and so despite the difference in free stream con- 
ditions between Young and Maas’s experiments and Hall’s theoretical 
problem, it does appear that the Young and Maas displacement is essentially 
a three-dimensional effect. An additional verification of this might be 
obtained from experimental displacement measurements using flat ‘two- 
dimensional’ pitots. 
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Finally, mention must be made of some recent experimental work by 


Davies (8), where the boundary layers on a flat plate and a cone in super- 
sonic flow are traversed by a pitot with circular cross-section. In all cases 
the stagnation streamline has a comparatively small deflexion towards a 
region of lower velocity. So far no theoretical verification or otherwise of 


this result has been obtained. 
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TWO-DIMENSIONAL CONTRACTING DUCT FLOW 


By J. C. GIBBINGS# and J. R. DIXON{ 


[Received 12 September 1955. Revise received 31 May 1956] 


SUMMARY 

This paper deals with the incompressible potential flow through two-dimensional 
contracting channels of finite length. These channel flows are obtained by first 
specifying flow patterns in the logarithmic hodograph plane. It is shown that 
certain of these patterns can result in infinite values, at points on the channel wall, 
of both the velocity gradient and wall curvature. A method of avoiding these 
undesirable features is given which alters the contraction boundary so as to replace 
them by wall portions along which the velocity gradient is constant. A numerical 
example is given. 


1. Introduction 

DESPITE a general increase in the velocity of a fluid stream on passing 
through contracting ducts in a potential flow, it is possible for the velocity 
to decrease along portions of the boundary. These adverse gradients can, 
in a real flow, cause boundary layer separation. In wind-tunnels this can 
seriously disturb the working section flow (1) which follows a contraction. 

Many methods are available for designing two-dimensional and axially 
symmetric contractions that have no adverse velocity gradients in the 
potential incompressible flow through them. It seems that all such con- 
tractions must necessarily have boundaries that are infinite in length. 
Goldstein (2) appears to have first provided a method of designing con- 
tractions having finite length boundaries. In his paper and in papers by 
Whitehead, Wu, and Waters (3), and by Lighthill (4), it is shown that 
adverse velocity gradients exist all along the upstream and downstream 
parallel walls. Goldstein further showed that these gradients could become 
infinite at points where the parallel walls joined the contraction curves, 
and for his particular design method he showed how these infinite gradients 
could be removed. 

It can be shown, by a proof similar to that in section 2 of this paper, 
that these infinite adverse gradients occur at each end of the contraction 
described by Whitehead; Wu, and Waters. 

Discussion of the relative lengths of different contraction shapes is of 
limited value if only the boundary shape is considered. For the velocity 
distribution across each end of a contraction only becomes exactly uniform 
at infinity up- and down-stream. In many real cases a desirable degree 
of uniformity can be specified and then the stations at each end of the 

+ Royal Aircraft Establishment. t The Bristol AerDplane Company. 
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ntraction where these distributions occur can be found. Then the con- 


traction length that is of interest is the distance between these stations. 


Thus both the centre line and wall velocity distributions are of interest. Un- 


fortunately both Lighthill and Goldstein discuss only the wall distributions, 


s they were largely concerned with the length of boundary wall. 
The design method used in this paper avoids these defects. It differs 
from Goldstein’s, and so, the occurrence of infinite velocity gradients is 


reconsidered more generally for applications using the present method. 


2. Preliminary general analysis 
In the present design method a flow pattern is specified in the logarithmic 
\lograph plane and certain features of these patterns are now considered. 
Writing 
ee dw ; 
Q) log j log g—10, (1) 
f 


where 


= velocity potential, 

stream function, 

y = coordinates in the physical plane, 

, 9 = respectively scale and direction of the velocity vector in the physical 

plane 

] orad d. 

Conditions in the z-plane at points corresponding to points of zero 
velocity’ in the flow pattern in the logdw/dz-plane are now discussed. 
We have 

dw “T7 2 
dQ ms (°) 
where U’ and V are the ‘velocity components’ in the Q-plane. 


From equation (1) we have 


hence lg = 4q(dQ+dQ); (3) 


ind also 
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where s is distance along the boundary. Then from equations (3), (5), and) and th’ 


(6) the velocity gradient is either 
dq g(dQ | dQ gu z 
ds 2\ds as) U2+y2 \ 
—_ ' ; " : , 1 
Similarly from equations (4), (5), and (6) the radius of curvature p is 
; : 
] dé 4h = (3 
pds UV : 
At stagnation points in the Q-plane, U = V 0, and in this limiting 
vent U .¢ aU 1e2U .. dU 
3 => if . £0 or pom aua if ; 0. 
U24+yp2 d\ 2 dV? d\ 
and also 
é wa j Lev xx : 
J pan a dV ae. d v f dV _ 
U2+ Vy? dl 2 dU? dl 
Thus at all points in the z-plane that correspond to stagnation points in | 
the Q-plane, the wall velocity gradient is infinite except when dU /dV = () 
whilst the curvature is infinite except when dV /dU = 0. 
Nu 


3. Numerical solution in the -plane for special boundary con- 


ditions boun 
The flow in the Q-plane is governed by either of the Laplace equations | 


insta 
V*¢ YO or Vd = Y, form 
The solution to be presented was obtained numerically and, as a conse: | 
quence, it was found more convenient to solve the second of these two 
equations. This was done by using Relaxation methods which, for a flow - 
past specified boundaries, are well known (5). “ 
The following method was used to obtain boundaries in the Q-plane, | a, 
which corresponded to boundaries in the z-plane along which the velocity | (0 
gradient had specified values. (C; 
The two boundary conditions are and. 
% = constant, (9 (12) 
aq = constant. (10 (d 
ds ; aqu 
Writing the boundary slope in the Q-plane as Y’ we have = 
Yy’ = LU (11) t 


U : radix 
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an and thus using equation (7) the boundary condition of equation (10) becomes 
either oa @2 
| WJ 17 on ; (12) 
‘ o(—@) (1+ Y’*)(dq/ds) 
ub . Y’q? 
on ; | oS (13) 
) 1S c(log q (1+-} 2)(dq ds) 
-8 
iting t 
° 
ta 
a’ 
ts L 
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? Fic. | 
-on- Numerical differentiation formulae provide a relation between the 
j . Ob Cus . 
oundary coordinates and the boundary values of or ¥ . For 


o(—8@) = a(loggq) 


nstance, with the mesh system illustrated in Fig. 1 and using one of the 
formulae given in (6), 
mise: f Dal; , 
bon A) ae a (1+-2r)po+ (1+1r)*by—r7hs, J. (14) 
flow 
The method of relaxation employed was to 
- (a) Specify a boundary shape in the Q-plane. 
i (b) Solve V2us = 0, for the boundary condition (9). 
(c) Numerically differentiate the specified boundary curve to obtain Y’ 
and substitute this together with the specified value of dq/ds into equations 
Q 12) and (13) to give = and op DS i 
o(—@) o(log q) 
i (d) Compute a corrected boundary shape from equation (14), which is 
' quadratic in r. 
(e) Continue the process until successive boundary shapes coincide to 
the desired degree of accuracy. 
(1] Alternatively one can obtain a boundary shape in the physical plane having a specified 


lius of curvature. This ability will be used in a future paper on aerofoil design. 
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4. Numerical integration to obtain the physical plane 
Numerical differentiation of the values of 4 obtained from the relaxation 
solution, substituted into the Cauchy—Reimann equations 
Cd ous 
@(logq) — @(—8)’ 
od os 
o(—8) e(logq)’ 


followed by a numerical integration will give values of ¢. From equation (1), 


i 
dz = —e8 dw 


q 
| (cos Odd — sin @dis)+-i(sin Odd 4 cos 6 dif) }. 
q 
Thus, 
dx = —(cos@0dd — sin @ dis), (15 
q 
dy | ain Odé _ cos 6 dys). (16 
q 


Numerical integration is conveniently performed along lines for which 
either logq = constant or #@ = constant. But even when computed values 
of é and & are known at regular intervals of log g and —@, equations (15) and 
(16) are still unsuitable for numerical integration. A simple integration by 


parts makes them suitable. For instance, when logg = constant, 


dx — —[d(¢cos 86) d(s sin 6)—{¢ sin 6+-% cos 6}d(—@)], 


l 
q 


dy [d(¢ sin @)+d(%s cos 6) + {d cos 0—y sin O}d(—4)]. 


l 
q 
Or, when —# = constant, 


dx cosoa()—sinoa(*) {ys sin 6 4008 6} d(logg), 
q 9. q 


: ; l 
dy = sin aa(* - cos d(*) -_{dsin 0+ cos 6} —d(log q). 
q q Y 
Simplification when also either 4 = constant or d = constant is evident. 


5. A particular contraction example 

Consider first the contraction as sketched in Fig. 2 where AB and EF 
are straight lines parallel to the centre line A’ F’, A and F being the points 
at infinity; BC and DE are curved boundaries along which the velocity is 
specified as constant; and CD is a straight line along which the velocity is 
to increase monotonically from C to D. 
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aie B 
A z- plane. 
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With this specification the corresponding flow pattern in the Q-plane is 


readily sketched as illustrated in Fig. 3, there being a source at A and a 
sink of equal strength at F. 





“a Log - sJ\ - plane. 
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Flows through similar contractions have previously been obtained by 
Hughes (7) and Waters (8). 
In Hughes’s chosen pattern in the Q-plane, both C and D were at 
# = oo and the source and sink were situated at B and E respectively. 
Thus the contraction obtained was of infinite length and had a singularity 
onsisting of a spiral at the point corresponding to C and D. Waters is 
nderstood to have obtained the contraction shape of infinite length 
orresponding to the pattern of Fig. 3, except that the source and sink 
were at B and E£ respectively. 
From the preceding results it is now known that along AB and EF in 
Fig. 2 the adverse velocity gradients will be infinite in value at B and E£. 
\s a first step towards a solution that avoids these infinite adverse 
gradients, the flow within a boundary of the type of Fig. 3 was obtained 
the relaxation method.+ The rectangular boundary is convenient for 
The solution of this pattern, given in Appendix I, requires less labour of computation 
in did the relaxation method. However, a later, modified pattern could only be solved 


the relaxation method, and it was thus more convenient to solve by relaxation from 


the start. 
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this method, hence the choice of the Q-plane in which to specify a flow | the regi 
pattern. Relaxation is also made easier by making the velocity drop from | this typ 
A to B the same percentage of the velocity at A as the drop from E to F} Thee 
is of the velocity at F. Then AB = EF and the pattern has an added 
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For this numerical example, the value of —@ at the points C and D was } pound. 

arbitrarily chosen as 1-0,¢ and the contraction ratio was specified as | contra 

exp(0-95) = 2-58, in order to fit a particular application, and to simplify | estima 

the relaxation mesh. | at eit] 

The solution given in Appendix I enabled a graph to be drawn illustrating plotter 

the change in contraction length with change in value of the percentage 

velocity drops. This graph is given in Fig. 4, where the change in length ' 

is expressed as a fraction of the contraction width at the low speed end. The 

From inspection of this, and again to suit the relaxation mesh, velocity , the be 

drops in the ratio of exp(—0-05) 0-951 were chosen. ey 

In order that the ordinates of the contraction shape should give a curve ae 

of sufficient smoothness a fair degree of accuracy was required of the } sionif 

relaxation solution. This accuracy was conveniently obtained using the add e 

definition of the relaxation residual as given by Woods (9). Relaxation in ” a“ 

+ It was thought that the success mentioned in the last paragraph of section 6 might there s 


not be achieved if this maximum slope was much greater. 
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' 
flow | the region of point A was carried out according to Woods’ technique for 
rom | this type of singularity (10). 
oF | The contraction shape obtained is shown plotted in Fig. 5 and the corre- 
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sponding velocity distributions along both the wall and the centre line are 
given in Fig. 6. The four infinite velocity gradients, two favourable and 
> two adverse, can be seen. 

To obtain a contraction shape along which the adverse gradients are finite, 
the relaxation method of section 3 was used to modify the boundary at the 
pints B and F of Fig. 3. The modified portions then corresponded to 
portions in the z-plane along which the adverse gradient had a constant and 
specific value. 


The method of Appendix II provided an approximation to the modified 


wa boundary shape required and also enabled an estimate of the change in 
ontraction length resulting from this modification to be obtained. This 
lify | estimate is shown graphically in Fig. 7 in the form of the change in length 
it either end, AS, as a fraction of the contraction width h at that end, 
Ing plotted against the non-dimensional adverse gradient, 
ALE 
h dq 
oti . 
q ds 
nd a - 
; lhere is difficulty in specifying a value of the adverse gradient so that 
ity 
the boundary layer in a real flow will not separate. Batchelor and Shaw 
were unable (1) to compute successfully the boundary layer development 
in : y 1a) 
Lh long a contraction where separation was occurring. In order that the 
significance of the scale of gradients in Fig. 7 can be visualized there is 


idded an auxiliary scale of cone angles of diffusers having corresponding 


" There appears to be a slight 1 in this paper at the foot of p. 182. The expression 
pa} I ! 
- here should read 


R F(4A’¢d A’¢,) —a*k, — J mr. 
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of a value of h dq A 
aly = — (0-15 
q ds 


in the numerical example. 
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velocity gradients. Hence there was only slight justification for the choie 
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The resulting approximate shape of the modified boundary in the Q-plane 
is plotted in Fig. 8. It is of interest in that an infinite adverse gradient is 
avoided, whilst there is still a stagnation point at H. As shown in Appendix 
II the slope is infinite at H. On moving away from H it diminishes from 
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this value with extreme rapidity, and as a consequence the relaxation 
solution was unable to distinguish, in the similar shape that was obtained, 
tween the true shape and the more common form of stagnation point 
that exists at a point of discontinuity in the boundary slope. No difficulty 
was, however, experienced with the relaxation solution and the final 
joundary shape was obtained quite rapidly after three successive 
pproximations from the shape given by the method of Appendix II. 
Fig. 8 provides a comparison between these first and final curves. 

The changed velocity distribution at the high speed end of the con- 
traction is compared with the original one, having the infinite gradient, 
7 D 
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in Fig. 9. This figure also illustrates the increase in length at this end of 


AS 
= 0-165, 
h 
compared with the approximate value from Fig. 7 of 
AS 
—— = 0-2]. 
h 


The discrepancy is not surprising, for the velocity potential at G in Fig. § 
can change quite rapidly with small movements of G along A B. 
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Though the change in the velocity distribution is seen to be appreciable, 
the corresponding differences in the contraction boundary shape were found 
to be not large. This thus suggests that in the potential flow the velocity 





gradients at the ends of the contraction are particularly sensitive to small | 


changes in boundary shape. 


The total length of the modified contraction is 1-28 times the width at | 


the low speed end. 


6. Practical application 
Whitehead, Wu, and Waters have shown (3) that a contraction shape 


obtained by two-dimensional flow analysis can be used as the shape for | 


an axi-symmetric contraction. They showed that gradients that were 
favourable in two-dimensional flow became more so in the axi-symmetric 
flow whilst adverse gradients were lessened. Analysis of the flow in square 
contractions (11) has shown that adverse gradients, when they appear, do 
so firstly along the corners. A comparison has been made (12) of the 
gradients along the walls of a two-dimensional, an axi-symmetric, and the 
corner of a square contraction, all having the same boundary shape. This 
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showed that the gradients were favourable in all three cases, those at each 
end of the square contraction being slightly less so than the corresponding 
ones for the axi-symmetric contraction. 

Though these are specific examples and not general results it seems that 
, suitably designed two-dimensional contraction shape can be successfully 
ised for axi-symmetric and square contractions. 

The difficulty discussed by Batchelor and Shaw (1), in predicting the 
boundary layer development along a contraction where separation was 
present, has already been mentioned. Present practice is to allow for 
viscous effects by adjusting the boundary coordinates, as obtained from 
the potential flow solution, by an empirical boundary layer displacement 
thickness. It is suggested that this can be more readily done with the 
present method where all velocity gradients can be made finite and of any 
desired value than with one in which infinite gradients are obtained. 


7, Conclusions 

The present paper has shown that a contraction of finite length must be 
wefully designed so as to avoid, in the potential flow through it, the 
nfinite adverse wall velocity gradients that occur in Whitehead, Wu, and 
Waters’ design (3). The elimination of these gradients is considered in 
greater detail here than in Goldstein’s paper (2) and a method is given 
whereby the value of the gradient can be more readily controlled than can 
be done using his or Lighthill’s (4) method. 
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APPENDIX I 
Solution giving the flow in the rectangular boundary in the 2-plane 
The origin in the {-plane, as illustrated in Fig. 3, is shifted to the centre of the 
rectangle by putting 


L Q—d— hry m+in (say). 


The rectangular boundary in this L-plane, as illustrated in Fig. 10, is transformed 
a unit circle, as illustrated in Fig. 11 by (13) 
sn(4AL)dn(3AL) 


cctenesist Al) 
en($AL) . _— 
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where the quarter periods K(k) and iK’(k) of modulus k and k’ respectively of the} around 
Jacobian elliptic functions satisfy the relations letterec 


K(k) Kk) A 



































> at =, A? 
B y 2 ” 
, : 1—en(AL) 
From equation (Al) 2 ——— , } 
; ( ¢ 1+en(AL) 
| ee 
whence en(AL) 
1+¢ j 
Fron 
vn ub 7 
¢ > N C-plane ) 7744! 
L- plane Ther 
of |e Na ty 
f x 
o ™ 7 - | and is 
, ‘. 
BA F € 
| Anap 
Fic. 10 Fic. 11 The 
} is shov 
: se ; . along 1 
and on the circle | =e, (A3 The fc 
whence en(AL) itanw. (Ad arms ¢ 
Along BE, L oe. Hence Fron 
en(AL) = en(Am—iK’,k) where 
en(iK’—Am,k) (7) an 
| 
i ds(—Am, k) 
—— a Along 
is a li 
and substituting into equation (A4) gives | dgjds 
The 
w tan iF ds( dm. k)] (A5 noted 
1 dq 
. 7 (he’ q ds 
and articular at B, tan—! ; d 
und in particular WR i (") of the 
The flow pattern in the ¢-plane, consisting of a source at A and a sink at F, 
given by w log[f- eilm+n -log[f- ei(-), (A6 
‘i Wri 
and substituting from equation (A3) this leads to the velocity potential ¢ as I 
§ = gig eee. (a7) B “1 bel 
" l -cos(w-+7) : Alo 


Specification of 8 and y enables K, K’ and k, k’ and A to be obtained from equations 
(A2). Then successive use of equations (A5) and (A7) gives the velocity potential 
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f the} around the boundary in Fig. 3. Using suffixes to denote values at the correspondingly 
lettered points in Fig. 3, the contraction boundary length stotq) 1s 
| D 
(A? Stotal | (gc—Pp) + | | p+ (be- dp) 
IB ¢ 2 VE 


2 ] 
dy, d »)—- — 7 5 
(p B) a4 709? $c)+ qr (bz—$p) 


l 2 
\(do—¢z) + (dp—¢¢).- 
q dptT Wd 


dB if 
From equation (A6) it is readily shown that by putting f = 0 along AF, then 
7 along ABCDEF, thus showing that the contraction widths at either end are 


ae 27 q4 and 277/qy. By arbitrarily placing A at the origin, then q4 = 1-0. 
Then the change of the contraction length, in terms of its width, with variation 
Nn q4/Up- 8 indicated by AS l 
a (dc—p) 
on qB 
and is plotte d in Fig. 4. 
APPENDIX II 
{An approximate solution giving the modified boundary shape in the ()-plane 
The boundary shape in the Q-plane, when the corners B and E (Fig. 3) are modified, 
} isshown in Fig. 8. An approximation is now obtained for the boundary shape HG, 
. long which, at corresponding points in the z-plane, the velocity gradient is constant. 
' The flow at this corner is approximated to by the flow due to a source at A, the two 
A4 rms of the corner extending to infinity. 
‘rom eqi oO 2 r 
From equation (2) dQ fi Vv 
= 9 
dw Q@' Q 
where Q is the ‘resultant velocity’ in the Q-plane. Substituting from equations 
anc 5) 
cade dQ 1dq (11 
U . 
dw q* ds qp 
\long GH, q will vary by only a few per cent. and so if this portion of the boundary 
sa line on which U/Q? constant, it will be approximately a line on which 
lq/ds constant. 
The dQ./dw-plane is as shown in Fig. 12, the points H and J being at infinity. It is 
A noted in passing, that point H corresponds to a stagnation point in the Q-plane, yet 
l de 
) 7, Is finite there. Writing p — dQ/dw, the p-plane is transformed to the upper half 
of the ¢-plane (Fig. 13) by the Schwarz—Christoffel transformation, which gives 
F, ee 
—— Ot-4(t+a)-(t+b). 
Ii 
A 
Writing ,/(6/a) = ¢ and substituting r = vt this integrates to 
Pp C[2r iva(1—c?){log(r—iva)—log(r+iva)}]4 C1, (Bl) 
\7 (, being the constant of integration. 
Along the line from G to J, 0 <r oo, and it is found that 
tions § " 


ss Na ; 
p C| 2r+2va(1—c?)tan“ +C. 
s 
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Thus C and C, are real, and hence at G 
aii? 
P C : —mva(c?—1) (B2 
and also at A where p = 0 and r = 7, (say) 
2ro+2va(1—c*)tan “+0, = 0. (B3) | 
0 
; d ( d ») - 
H H, L du fog <* = p- plane. | 
Vv 
* 
} 
m | 
oj/A I 
G ° 
Wee 
Fic. 12 
t - pl ane. 
I MH G A I 
wee. 
a 
~~ ! 


Fic. 13 


Along the line from H to M, iva < r < ivb and writing r ir’ 
—C isd r’—Na 
f —- 2ir’ + iva(1—c?)log ———. (B4) 
( r’+~wva 
- p_. nO, seo, 
Thus C, 0 and at M iva} 2c+(1—c?)log mit 
C+ 


Comparing this equation with equation (B2) shows that the ratio of the values of ; 
at G and M is independent of a and C, depending only upon c. As only this ratio 
and the value of p at G will be varied, the value of a can arbitrarily be put equal 
to 1-0. Thus equation (B1) becomes 
) P , r—t 2 
P 2r + 7(1—c?)log —.. (Bd 
ys 


The flow in the ¢-plane is that given by a source at A, thus 


w = log(r?—7?); 


dr r2— 2 - 

hence = = ; (B6 
d ( =) _nr—rid (1 dw 
anc Pp iio eg a = a , og - 
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| Substitution into equation (B5) gives 


l d du 2r2 r r—1 
log +4 (1—c?)log —.. B7 
B2 2C an dz r2— re ‘72 ra ii r+ ~~ 


integral of the first term on the right-hand side is 





: r—To 
al, 476 log (B8) 
B: 7m 

integration of the second term leads to an integral of the type 
, log(1 7) ° 
: [ way, 
, J 
0 
that is, Euler’s dilogarithmic integral. So equation (B7) is integrated numerically. 
{long the boundaries to be considered r is either wholly real, when 
7 i : l 
log 2itan-!-, (B9) 
7 i rT 
s wholly imaginary, when 
? a , l—r’ 
| 1g . im + log - (B10) 
l a l+r 
a : 
alue of (B8) becomes infinite at A, thus in integrating from G to A, the whole 
f equation (B7) is integrated numerically; thus, noting equation (B9), 
ro 
} l r a ] 
logqg,—log qd, 5 ~| 2r + 2(1—c?)tan! dr. (B11) 
2( 1 k r2 r2 r 
0 
limit of the integrand at A as r > 74 is 
ro ] 
I 2 aad ee WP as 
} 1 +75 tan—1(1/ro) 
Substituting equation (B10) into equation (B7) gives the ordinates of the curve GH as 
| dw ¥ 7 { 1—r’ | 
log | .| 2r+-7(1—c?) log --Ligr} | dr. 
(B4 sal dz J, | ye — | ey} - } 
reduc 
I ’ ‘ ios » 
= log q hor c?)log(r’2 + r2)+ constant, (B12) 
y 
| c il l—r ‘ . 
' 6 2\? r, tai | +(1—c?) - ,, log dr’+constant. (B13) 
q / e t’ | ret r’2 1+’ 
0 
Specifying values of 7,, equation (B3) gives corresponding values of c. Use of 
tions (B11) and (B12) gives (1/C)[logqg,—logqy], then specifying logqgy and 
ting q4 1-0 gives C. In the present case ob 0 on the contraction centre line, 
{7 on one wall, so that the contraction width at the low speed end is 27 and equa- 
h de . , 
' n(B2) then gives { t) . Values obtained when gy 0-95 were 
q ds/, 
(Be 
, 0-] 0-2 0-5 1-0 2-0 


- ] — —()-16126 0-08258 —0-1439 —0-2548 —0Q-4912 
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A difficulty arises in computing the coordinates of the curve GH, for the integrand 
in equation (B13) tends to infinity as 7’ > 1-0 at H. This is avoided by computing 
from the following series expansion the integral 








e , 


r , , 
5 “ali—r’) or 
ror 


all l—rgl—r’, 1—3rf (1-—7’)? . } 
+All 14+R 4 TR? 9 UT | 
{, ,1l—rgl—r’  1—3r3 (1—7r’)? \| 

le a—r)| 1+ _— -1- ~ wt (8 

Bt MU tite 2 tase a thy OM 
As r’ +.1-0 the right-hand side tends to zero. 

To investigate in more detail the shape of the curve GH near to H, equations | 
(B12) and (B13) are, noting equation (B14), expanded in the series form 








(1 ] ) kar(1—c? | . 1—r’) sis 
; (log q ogq $77\ 1 — C* Sera om 7 T eeely 
5q (loga — logqu) = 4 Nia! Gabe 
Writing 
l (1—c?)(1+log2)—2 ; 
g0(— 0+) = [r+] + 
+log(1—r’) = ae ie ed Values 
0 
Including only the terms in (1—7’), 
6+-6 ] 2 
mo td Ss [: + log 2— = +log(1- r)). 
logg—logqy 7 1—c2 1.G 
This illustrates the extreme rapidity with which the slope of the curve changes fron 2 me 
its value of — 00 as a point moves away from H. aa 
‘a " ; - re 
The effect, on contraction length, of the modified boundary 
To obtain an approximate value for the change in contraction length due t 4M 
rounding the corner B in the logarithmic hodograph plane consider the two flow re 
boundaries of Fig. 14, that of case II being the present approximate one. | 5s R. 
In case II, on going from H towards J, the influence of the rounded corner GH 6 C. 
will have a lessening effect and the flow will tend to that in case I. : | 
Thus the difference in contraction boundary length between case I and case II : 
will be a tn 7H 
; . - l 
AS = | 7 $— | (b= [¢n,—$e,,| 
: ‘ B 
7 / 8. M 
In case I the flow in the Q = log(dw/dz)-plane is given by | 9. | 
w log(Q—a) — log(Q+a). | 
Q—a 
Hence e” = . 10. 
Q+a 
ll. J 
At B where Q = 0, « Lo = 4, 
12. 
dp, 0 
In case II, from equation (B6), 
13, | 
elt 72 rz 


At G where r = 0, e" —riu = 7, 


and hence $6) 2log 7. 
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I 
I 
4 | -8 
H 
| 
I B 
A GA lo 
- oe a an log + 3 % 
CASE I CASE IL 
Fic. 14 
P : 2 
Putting gy = 9-95, AS = 995 °S"e 
Writing h contraction width at low speed end, 
AS log rs 
h 0-957" 
Oy aes 
Values obtained are plotted against - in Fig. 7. 
ds 
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A NEW EXACT SOLUTION OF THE EQUATIONS OF 
VISCOUS MOTION WITH AXIAL SYMMETRY 
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[Received 16 December 1955] 


SUMMARY 


An exact solution for the motion of viscous fluid flow for axially symmetric motion 
has been found by Squire (1) in the form % = vrf(y), where (7,6,¢) are spherical 
polar coordinates, pu cos 6, and v is the kinematic viscosity. Squire has shown that 
a special case of his solution gives the flow for a jet emerging from a hole in an infinit. 
plane wall. In the present paper, another exact solution of the Navier—Stokes 
equations is obtained for the steady flow of a viscous incompressible fluid for th 


case of axial symmetry in the form & = vr"f(8) provided that n 1, 2, or 4, 
For n 1 the solution becomes that of Squire. For n 2 we obtain the solution 
% = vr?C,sin?@ (C, being a constant) which is well known. The case n = 4 appears 


to be new. 


1. The equations of motion 
WE take spherical polar coordinates (7, 4,4) with @ measured from the axis 
of symmetry. The velocity components are (wu, v,0) measured in the direc- 
tions (r, 8, d) respectively, and all the quantities are independent of ¢. The 
equation of continuity is (2, § 41) 

ap 7) + ~ie 0 piersin es, 


and the equations of motion are 


cu veu wv 1 Op > 2u 2év 2vcoté 
e¢—+—-—_—. —-— —+y V¥u———.— = — - : 
cr roe r p or - Fe a 
ov vov. w 1 op ; 2 eu v 
u 7 + i i le Doel i +v V70+ 9n0 .9..:..90)? 
or rob r proe r 06 r* sin?6 


where p is the pressure, p is the density, v is the kinematic viscosity, and 


l ¢ 0 l e 0 
V3 2—} 41 —___ _ |gin 9 _}. 
r? Or ( *| r? sin @ 00 (: a 


We shall assume that % (the Stokes’s stream function) has the form 


ob — vr"f (8), (1) 
so that He 1 eb prn-2 
~  - rgsin@ red sin 6° 

1 edb m rnr—2 
v= y = f (8) 


rsin 6 é@r sin@ ° 


f’(9) 


(2) 


(Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 1 (1957)] 
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The equation for the stream function ¢ is (2) 


l c (ub, Ds) 2D s a0 sin 6 Oys| 


vDAd, (3 
r=sin@ a(r,@) r?sin®0|— r 06) ¢ 
os sin8@ ¢ 1 ¢ 
here D? pe r . 
— or? a land a} 
Substituting the assumed form (1) for x, 
D*y = vr"-*( f”— cot 0f’+n(n—1)f ); (4) 


and hence from (3) we get, after some simplification, 
r®-1[ _ 2n?(n— 1)cot 0f?+ nf 4n—3-4 3 cot?6} ff’ 
3n cot Off" + (n—4)cot 6f’*—(n—4) f'f"+nff” 
sin 6[n(n— 1)(n—2)(n—3)f—(2n?—6n+-9+ —_ cot 6 f’+ 
(2n?—6n+8-+3 cot?@) f”—2cot6f”+f'*], (5) 
where dashes denote differentiation with respect to @. 


2. Equation (5) is the differential equation which f(@) must satisfy in 
order that vr”f(@) may be a possible form for the Stokes’s stream-function. 
Putting | in equation (5), we can show that it has a solution 

2 sin?6 
C+1—cos0° 
The corresponding stream-function is & vr f(@). This is the case investi- 
gated by Squire (1, 3) 

If n | the presence of the term r”-! in (5) shows that the two sides of 
this equation must separately vanish. We shall find a solution which is 
common to the two differential equations thus obtained. 

We find by actual substitution and simplification that the left-hand side 
of (5) equated to zero is satisfied by f(@) = C,sin"@, where C, is a constant. 
However, when we substitute this value of f(@) in the right-hand side of 
5), the latter reduces to 

n(n—2)*(n—4)C > Sin" 36, 

It appears, therefore, that n must be 0, 2, or 4 in order that C,sin"@ may 
bea possible value of f(@). Asn 0 is trivial and n = 2 gives a well-known 
solution, we consider ‘ied only the case n = 4. Putting n = 4 in (5) the 
equation becomes 

+ f[ ff” —3 cot Of” + (13+ 3 cot?6) f’ —24 cot Of ] 

sin 6| f’*Y —2 cot 6f” + (16+ 3 cot?@) f” —(17+-3 cot?@)cot 6f’+ 24f] 


2 ° d r 1 9 “7 
| 47°f—sin 8 797 cot 6] tu 3 cot 6f” +-(13+-3 cot?@) f’—24 fcot é] = 0. 
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Hence any solution of 
f” —3f" cot 0+ (13+-3 cot?@) f’—24fcot @ = 0 
gives a possible stream function. Now 
f” —3f" cot 0+ (13+ 3 cot?6) f’—24f cot é 


d ‘ ” ’ 9 
(j4—2Zeot Ns —f’ cot 0-+12f). 


Hence we require a solution of 
f’"—cot 0f’+12f = Asin*6, 
where A is a constant. Now the general solution? is 
» a , ' , 
f (0) = sin®6[},A + BP4(n)-+ CQ5(u)], 
where P3, Q3 are the derivatives of the Legendre functions P, and Q,, 
p = cos@ and B and C are constants. 

A particular case of the above solution is obtained by putting C = ( 

and A/B = 15; we thus get 
w% = kr*sin?@ cos?0, 
where k is a constant. 

This solution is of interest since it makes u = v = 0 for 6 = 47, and this 
represents a motion with a plane boundary. The stream lines are rectangu- 
lar hyperbolas in this case. 

This may be regarded as the motion in a slightly divergent jet striking 
a distant wall at right angles, but as in the case considered by Squire, the 
boundary conditions are not satisfied at the outer surface of the jet. 

It is easy to show that in this case, at any point (r, 0,4), the components 
of vorticity are . 

y f= _= 6, C = 2krsin 8, 
and q? = 4k*r* cos?6, 
where q is the resultant velocity. 

In conclusion I thank Dr. Gorakh Prasad, Allahabad University, Allaha- 
bad, under whose guidance this investigation was carried out. 
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THE LAUNCHING AND PROPAGATION OF ELASTIC 
WAVES IN PLATES 

By H. PURSEY (National Physical Laboratory, Teddington, Middlesex) 

Received 16 February 1956] 


SUMMARY 

Integral representations are derived for the displacement at an arbitrary point 
n an isotropic plate due to a stress distribution on the free surfaces which varies 
sinusoidally with time. It is shown that evaluation of the integrals leads to secular 
equations identical with those originally obtained by Lamb, and the solutions of 
these equations have been investigated for the symmetric and antisymmetric wave 
types in plates up to 8/7 compressional wavelengths thick. New tables of solutions 
are provided for the two principal waves and graphical representations of the solu- 
tions corresponding to complementary waves are given in Figs. 1-6. Difficulties 
which arise in the use of Lamb’s tables have been avoided by the method of normali- 
zation adopted in the present paper. 

Formulae are derived from which the amplitude and group velocity of a given wave 
may be calculated, and their application is illustrated by means of a worked example. 


1. Introduction 

THE first rigorous analysis of wave propagation in an isotropic plate was 
given by Lamb (1) for the two-dimensional problem of a harmonic wave 
travelling in a direction parallel to the medial plane. He considered two 
types of wave, one symmetric and the other antisymmetric with respect to 


he medial plane, and derived equations relating their phase velocities to 
the thickness of the plate. It is apparent from these equations that the 
phase velocities are many-branched functions of plate thickness, and with 
the exception of one branch in each case their values are complex when the 
plate thickness is small compared with the length of a plane transverse wave 
in the material. Hence, only waves corresponding to one symmetric and 
one antisymmetric branch can propagate in a thin plate. These are referred 
to as ‘principal waves’, and corresponding solutions of the secular equation 
are tabulated by Lamb. 

The symmetric modes of propagation have been studied by Holden (2), 
who shows how the approximate form of the solutions may be deduced from 
Lamb’s secular equation without undertaking a detailed numerical analysis. 
Holden also considers the case of symmetric waves propagated in cylinders, 
while the principal antisymmetric branch has been studied by Osborne and 
Hart (3) for plates, and by Hudson (4) for cylinders. 

The purpose of the present paper is to show how the amplitude of the 
listurbance is related to a given stress distribution on the free surfaces of 
the plate, varying harmonically with time. Two types of source are con- 
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sidered, one giving a two-dimensional field, as in Lamb’s paper, and the 
other having circular symmetry about an axis normal to the surface of the 
plate. The method of solution follows that used in the analysis of the field 
from a radiator on the free surface of a semi-infinite solid by Miller and 
Pursey (5),+ leading to definite-integral representations of the field com- 
ponents at an arbitrary point in the plate. The value of the integrals is 
shown to be the sum of the residues at the poles of the integrand lying in 
the upper half-plane or on the negative portion of the real axis, and investi- 
gation of the positions of the poles leads to the secular equations originally 
obtained by Lamb. New tables are provided, showing the velocity of 
propagation as a function of plate thickness for the two principal modes 
in materials with Poisson’s ratios of } and 4. The parameters are normalized 
by comparison with the velocity and wavelength of a plane compressional 
wave in the material. This avoids difficulties which arise in the use of 
Lamb’s tables, where the plate thickness is measured in terms of the length 
of the appropriate principal wave—a quantity which is itself a function of 
plate thickness. 

Consideration is given to the complementary modes of propagation (that 
is, modes other than principal modes) in both the symmetric and anti- 
symmetric cases, and Figs. 1-6 show the relationship between phase velocity 
and plate thickness for propagating waves in plates whose thickness does 
not exceed three compressional wavelengths. 

Formulae are derived from which the amplitudes and group velocities 
of the normal surface components of all propagating waves may be obtained, 
and their application is illustrated in a worked example in which the plate 
is assumed to have a Poisson’s ratio of 4 and a thickness equal to 2/7 
compressional wavelengths. 


2. Derivation of integral representations for the displacement 

components 

It is convenient to reproduce here the equations of displacement, stress, 
and propagation derived in section 2 of Paper I, which also find application 
in the problems considered in the present paper. Cartesian and cylindrical 
polar coordinate systems are used, the z-axis being in both cases normal to 
the medial plane, which also contains the coordinate origin. The relevant 
equations are 


tp = — 28 (eOE iM), tap = — (TE iets), (0 
pw” dz pw*\ dz 
2 2W, 
ih ¢?—ki)A, = 0 ; Wr — (C?—k3) Wy = 0, (2) 
dz? dz” 


+ Subsequently referred to as ‘Paper I’. 
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ru 2 _ d*A , ° LW; 2/..2 2 , 
Oe Bp = — pF + 2 F + pu? 2)0Ay, (3) 
C44 - - 
ee Say = — (Me + air Se + omy) (4) 
ci, dz? dz 
» yee Cy (AWa, P “ 
ee as (u? = Wis; Upp, = — | = one tn) (5) 
Awe ((2_13)Ap,—0, “Nes_ (2 1)Wy, = 0, 6) 
dz? dz* 
~ Ep — 4 d*Ap, L o¢ tn, eS p?(p?- 2)C?Ap, 9 (7) 
Che 70 : dz? dz ¥ 
Ww. (d*W,, l yoy 
= Fy, — ure An +o), (8) 
C44 az dz 


where 
u,, U,, and wu, are components of displacement, 
2, 2, and 22 are components of stress (Pearson’s notation), 
p is the density of the plate, 
w is the pulsatance of the vibration, 
c,, and ¢c,, are respectively the compressional and shear elastic constants 
of the plate, 


and ky = wy(p/Cy), ky = wal(p/C44), (9) 
A=V.U, W VxU|, (10) 
iv / (Cy1/Ca4) ky k. VJ{2(1—o)/(1 —2o)}, (11) 


ais the Poisson’s ratio of the plate. 

The suffix F, By, or B, indicates that a variable has been transformed 
by the Fourier or Hankel integral transform, defined respectively by the 
equations 


x 


QP} (Z) [ d(ax)e—*s* da, pz,,(S) [ d(r)rJ,,(Cr) dr. (12) 
“a 0 


As in Paper I, it will be convenient to normalize the equations by 
measuring lengths in units of 1/k,. Equations (2) and (6) then become 


o : ¢° JA; U, —— (¢? pw) Wp 0, (13) 
dz dz? 

ew (S° 1) \; 0, — (¢? pw) We, 0. (14) 
dz2 dz? 


the remainder of equations (1) to (8) being unchanged, provided that all 
symbols are now regarded as denoting quantities measured in normalized 
units. 
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Case (a) Field invariant with y 
SYMMETRIC WAVES 
It is convenient to introduce here the notation 
S” = sinh{z,/(¢?—m?)}, Cr = cosh{z,/(¢?—m?)}. (15) 
The appropriate solutions of equations (13) for waves having symmetry 
with respect to the medial plane are then 


Ar = A, Cc}, Wy = B, SH, (16) 


Substitution in equations (3) and (4) leads to the following expressions for 
the transformed components of stress: 


w ¢ 9 9 - 

ae se A, p?(20?—*)¢ 2 20B, C,/( (f?—p?)Ce (17) 
C44 

a “yp 21 A, ply (C?—1)S7+ B,(20?—p?) SE. (18) 
44 


If we now write for the transforms of the symmetric components of 


stress on the free surface 


Ail) (22,7), a? g,(S) = (2%). a’ (19) 


where a is the half-thickness of the plate, we obtain for A, and B, 


ie eee 27 ) SH. —p2)g,(£)C#) 
aan ptc?, nian" we) f(C) Sh + 20ly/(L?@—n*)g, (2) CH} 
By =P" _ (272 pg, (0)C1—2it (C21) f(C) SY, 


~ 2, D(Lsa, p) 
where 
D(C3a,p) = (202@—p2)2S# C1402) (22—p2),/(Z2—1)SLCH. (20) 


Using equations (1) and (16), we then find for the transforms of the field 
components 


l 
l ‘ i es f(2¢2— u?) S# S1— 22281 Se) + 
Us} Sus Ditsa,n) A(S){(2 > be » C Z C Zs 
tiLgy (£){2(L2—p2)/(2— 1) Ce S3—(202—p2)CL SH], (21) 
1 *y ¥\C 9 > Y 2 2 / 2 
Mee = > — [ib (0 (20°?) Se C}—2y(22—p2) (21) 83 C8} — 


—(l?—p*)g,(f){202C# CL—(20?—p2)CLC#}]. (22) 
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ANTISYMMETRIC WAVES 


The appropriate solutions of equations (13) for waves antisymmetric 


with respect to the medial plane are 


A,y=A,S}, Wr = B,C. (23) 


If we write for the transforms of the antisymmetric components of stress 
the free surface 


Fa(S) C2 )2— a> Jo(C) (22). 


2 (24) 


a? 


nd proceed as in the case of symmetric waves, we obtain for the transforms 


f the field components 


ata V(¢ 1) fo(f){ (202 — wp?) C# C1—207C1 Ce} 
Cag SS, % b 
UCgo(C (24) (C° 1°) ( (¢ 1) S# C} (2¢?- -p?)S7, C#}). (25) 
| ‘ 72 2\(% VI » v2 2 v2 C1 Sei 
} D : WJ o(¢ \<$" pe") a”™2z Hy (C L )l(G 1) a'zs 
O44 (S34, | 
(C2? — p*)go(C){ 207, S# S} (2¢?—p?) 2) $1 S#'], (26) 
ynere 
Dy: a, pu) = (202@—p2)2S81 CH#—402,)(02—p2),(C2—1)S# CL. (27) 


h) A rially symmetru field 
SYMMETRIC WAVES 
The solutions of equations (14) corresponding to waves having sym- 
etry with respect to the medial plane are 
Ap A, C1, W,, = B,S#. (28) 
If we write for the transforms of the symmetric components of stress on 
the free surface 


f(C) C2 BRB Jz=a> 93(C) (27 ,,)- a (29) 


d proceed as before, we obtain for the transforms of the field components 
! /(C2—1) fa(C){(20? — wp?) Sh S!— 202,82 Se} 


Cga(C){2y (C?—p?)y (C2— 1) Ch SE — (20?— pw?) C7 SB}], (30) 


By = pee IS falM2(O—HP dy (2 1) Si Cl — (20 p? Se CB + 
; (0:4, 2) 


L*)gg(C)i (20 p*)¢ a ( . 
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ANTISYMMETRIC WAVES 
The solutions of equations (14) corresponding to waves antisymmetric 
with respect to the medial plane are 
Ap = A, 83, Wy, = B,C. (32) 
Writing for the transforms of the antisymmetric components of stress on 
the free surface 
Fa(S) (22 p,). ao g4(C) (27 p,)- a (33) 


and proceeding as before, we obtain for the transforms of the field com- 





ponents 
] 
bes = = Cyt (2¢2—?)C# C! 2¢2C1 Cr} 4 
=Bo owe D,(C;a I Lv > hs (¢) > I ) “sz > zs 
+ Lay(C){2,/(C2—p?),/ (C2? — 1) S# C1— (20? 2) S1 C#], (34) 
l 


UrB, or [ofal SHA (CP — pw?) (CP— 1) CG SE— (20?— pw?) CE S+ 
C44 D,(C: 4, 1) 
+ a/(f?— pe? )gg(C){(20? —?)S3 S#— 207.8 ST]. (35) 
The displacement components are obtained by application of the appro- 
priate inverse transform formulae. 


C'ase (a) 
SYMMETRIC WAVES 


cO 
l r 9 J a ¥ “9 » ’ ’ ye ’ Y 
w= 5—— | IWC 1A M{22—p2) Se St 20282 Se} 
“71C 44 
— Re aes ee "ee a eibr dl - 
+ 409, (C){24/(C?—w*),/(C?— 1) CR S}— (2? x?) C2 SHY] (36) 
D\C;a,p)’ 
l q *Y J Y { » “9 9 Yu 1 » ‘oO » 9” 1 ya) 
u, = — [eC fy (C)} (20? — px) SH ¢ = 2,/(6° bm )y (CP 1)S2 ¢ ef 
=7C 44 , 
¥9 9 Y\(o 26m 1 972 1 1) ese dl o= 
V(C°—p*)gi (Ch 207C# Cl —(20?@—p?)C} CH} lj -. (37) 
W(C; a, 1) 
ANTISYMMETRIC WAVES 
I 1 r2 ] -(7\S (972 2 (#1 9721 (CH 
u 2mC44 LV (¢ falC)i(2¢ b ) a 2S a on 
ve 
*y y ¥9 9 v9 ’ ” ¥9 9 ’ ’ etst dl ° 
+l go(C){2,/(C?—p?),/ (C2? — 1) S# Cl—(2¢?—p?)S1 C#}] ———, (vd) 
; “D(C; a, 1) 
I 1 rt (7\S (972 2\(H% VI » v2 2 v2 1 Op 
Us ? ae [aC fa(C){(2¢ B ) a”z mn iG" [es (¢* 1)¢ a St f 
=7C 44 
v9 9 Y\(o 72 Ou 1] v2 1 4) else dC 2g 
/(S* BM )Jo( Si 2C°S# a; (26° 2) SI SHY] (ve 


D(C; a, yu) 
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(ase | b) 


SYMMETRIC WAVES 
u, = (C21) fa(){(242 pp?) St S1— 222.8] Sv} 


Cag(C){2./(C2—p?),/ (C2 — 1) C# S1—(202?—p?) Ch S#}]— » (40) 





[C fa(){2,(C2—p2),/(C2— 1) 8} C#— (202—p2)S# CH 


tw 





(f2— 2) 22 p2)C1 Ce— 2008 C1) CI, (or) dé 
S [4 j 2 “> ( 2 &S C zs Pr ° 
. site, ; D,(f;.a, 1) 
ANTISYMMETRIC WAVES 


ul (C2—1) f,(C){(20? pp?) C# C1— 202C1 C#}4 


a 


, eer 
(q,(C)§2./(C* 1“), (C* ] Se C2 (2¢? 12)S1 ( hi} = sh eos. (42) 
1 VS PUINNS 2 I t WD Ga,p) 


(21) C1 Se — (22—p*)C# $3} 


r 
- 
to 
2 
> 
t 
I~ 


- — ee we ces cnra aler) & 
V(S2?—p?)ga(C){(202?— pw?) So St — 202.8 SI}] . (43) 
D,(C; a, 2) 
[he integrands in the above expressions can be rearranged so as to 


nvolve only even functions of the radicals ,(¢?—?) and ,(¢?—1), and 


therefore have no branch-points. However, zeros of D,(¢; a, ~) and D,(¢; a, x) 


give rise to poles on the real axis, and consideration of the effect of finite 
lissipation in the medium on the positions of the poles leads us to introduce 
ndentations of the contour, as described in Paper I, section 3, passing above 
the positive and below the negative poles. 

lhe integrals corresponding to case (a) are convergent at infinity in the 
pper half-plane, and are therefore readily evaluated by the theory of 
residues. To evaluate the semi-infinite integrals which arise in case (6), we 


nake use of the following results from the theory of Bessel functions: 


J, (2) fH) (z)— HW) (zexpiz)}, 
J,(z) fH) (z)+ HY (zexpiz)}, (44) 


where H\) and H{ are Hankel functions of the first kind (see 6, ch. 3). 
lhe integrals have one of the following two forms: 
T= | xilSo(Sr) dl, Ty = | x00) (Er) aE, 


0 0 
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where x, and x, are respectively odd and even functions of ¢. Applying the 
formulae (44) above, and changing to a new variable of integration, —:, 
in the second term of each integrand, we readily obtain the results 

x a 


I, [ xi(S)AG? (Er) dé, I, 5 





x2(C) M4? (Sr) dC, (45 


L x 
the contour passing above the branch-point at the origin in each case, 
The integrals may now be evaluated by the theory of residues, as in case (a). 


3. Determination of velocities of the principal waves 
The affixes of the poles of the integrands in expressions (36) to (43) are 
determined by the equations 
D(C: a, p) 0 (n 1,2). (46 
We exclude the solutions € = p, n = 1, and [ = 1, n = 2, inspection of 
the integral formulae showing that these do not correspond to poles of the 
integrand. Since (46) is an equation in @?, the poles are distributed sym- 
metrically about the origin, but only the negative poles are enclosed by 
the contour of integration. If the positive poles, arranged in increasing 


order of magnitude, are denoted by p,,,,,,m 1, 2,..., then the corresponding 
phase velocities are given by 
_ am V. Pram (47 


€ 


Rearrangement of the expressions (20) and (27) for D,(¢: a, u) leads to the 


where V, is the velocity of a plane compressional wave in the medium. 


following secular equations: 
Symmetric waves 


(2p?—p*)? (p> — pw?) (p?— 1) Sa(p)Ch(p) 


> = =" - (48) 
4p° Sh(p)e 1 P) 
Antisymmetric waves 
(2p?—p?)? __V(p?—p?)./(p?— 1) Sh(p)CU(p) (49 
4p* Si(p)Ch(p) 
We next consider the asymptotic solutions of equations (48) and (49) as 
a—>cxandasa-—. Provided |f| > p, we have 
SiH seo} ; 
lim —*—? — lim —2—4 = ]. (50) 


a 0 SH C1 a-—« S1 Ce 
The limiting equation for both symmetric and antisymmetric waves 
5 ©&Y p A 
with |C| > p is therefore 
I2— 42)2 2 /(m2t— 2 2 51) 
(2p?—p?)? = 4p?,/(p?—p?)\/(p?—1), (51 
which is the secular equation for surface waves propagating in a semi-infinite 
medium. 








Whe 
ofthes 
from si 


long in 
hence 1 


For. 


Hence. 


there | 


Equati 

32, an 
equal t 
include 
in this 


4, Pre 

It is 
the ra 
antisy! 
0< (C 
imagin 
is clea 
possib! 
compr 
comple 
A. N. 
menta 
each b 
that t! 
that t 
in the 
would 
attem 
preser 
findin; 
locati 




















ON ELASTIC WAVES IN PLATES 53 


When |¢ u the limiting formulae (50) are no longer valid, and solutions 
of the secular equations exist which correspond physically to waves reflected 
from side to side across the plate. The path lengths of such waves will be 
long in comparison with the distance travelled along the plate axis, and 
hence they will in practice be attenuated rapidly. 


For small a we may make the following approximations: 


(205m—H?)? = 4Pin(Pim— 1), 
4 344 
; U Sy 
Hence, 4 3 4 9 
Pil 4(u°— 1) P21 4a? (u* 1) 


there being only one real solution in each case. Therefore, for small a, 


= 4( pu" Laie ‘ : 4a?(u? 1) . bs 
y2 Ls oe ees (52) 
” u4 = 34 


i 


Equations (46) have been solved numerically for values of a between 0 and 


» 


2, and for p v3 and 2, corresponding to plates having Poisson ratios 
equal to | and 4 respectively. The results are presented in Table 1, which 
includes an auxiliary table of p,,va for 0 < a < 1-2, to assist interpolation 


in this range 


4. Propagation of complementary waves 

[t is readily established that equations (46) have only one solution in 
the range p C ©, corresponding to the principal symmetric and 
antisymmetric waves discussed in section 3. However, in the range 
0< \C pu the arguments of some or all of the hyperbolic functions are 
maginary quantities, and the behaviour of the functions is oscillatory. It 
is clear from the limiting expressions in section 3 that no solutions are 
possible in this range when a is small, but as a approaches the length of a 
compressional wave in the plate, solutions arise which correspond to 
complementary waves, characterized by high phase velocity and dispersion. 
\. N. Holden (2) has investigated the behaviour of symmetric comple- 
mentary waves in considerable detail, and has found confining lines for 
each branch of the solution of the secular equation, thereby demonstrating 
that the branches do not intersect. It is also evident from Holden’s work 
that the curves representing each branch are very irregular, particularly 
in the region near p |. To provide an interpolable table of solutions 
would therefore entail a large amount of computation, and this is not 
ittempted either by Holden or, except for the principal waves, in the 
present paper. Instead, Holden has studied the symmetric branches by 
finding a set of smoother curves about which the branches oscillate, and 
locating the points where the branches intersect these curves. In the present 
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ON ELASTIC WAVES IN PLATES 
[It is convenient to express the secular equations in the following forms: 
0) p | Symmetric waves: 
tan a,/(u?— p?) 4p*, (u?— p*),/(1—p*) (53) 
tana, (1—p? (2p?—p?)? 
Antisymmetru waves 
tana,/(1 p=) 4p, (pu? ] ? - 
A Ll fs ae t= lat tJ (54) 
tan a, (u-— p* (2p*— *)* 
Pp ri Symmetric waves 
tan a, (u"— p*) 4p, (ue py (p? 1) (55) 
tanh a, (p?— 1) (2p?—p*)? 
Antisyn metric waves 
tanh a, (p?—1) 4p,/(u? — p*), (p*— 1) (56) 
tan a,/(u?— p*) (2p?—p*)? 


For the numerical solution of equations (53) and (54) we make the change 


of variables 


» 9 9 9 
’ a, (1—p*), m (uw? —p*)/./(1—p*) 
und express the right-hand sides as functions of m and yp, so that, for 
example, (53) becomes 
tan mx . — 
F(m, p), (57) 
tan x 


where F is independent of x. For any given m we can calculate F and hence 
solve (57) for x. In particular, if m is an integer tan ma/tanx reduces to a 
polynomial in tan a. 

Similarly, when | < p < p» we may first transform equations (55) and 
56) by changing to variables x’ = a,/(p?—1) and m’ = ,/(u?—p?)/,/(p?—1), 
ind solve for 2’ in terms of m’. When 2’ is large a good first approximation 
s obtained by replacing tanh 2’ by unity. 

Figs. 1—6 illustrate the symmetric and antisymmetric branches for values 


of a between 0 and 8, and for pu V3, 2, and v5. 


5. Calculation of the group velocity 

The rapid variation of p with a indicates that the plate is a highly dis- 
persive medium, and therefore in order to calculate the time of travel of 
\ pulse of energy between two points on the surface we need to know the 
‘group velocity’, which is the velocity at which the pulse envelope travels, 


is distinct from the phase velocity calculated above. A detailed analysis 


a] of propagation in a dispersive medium has been given by Sommerfeld (7) 
bl ind Brillouin (8), and is summarized by Stratton (9). A rigorous treatment 


if a given problem would involve a formidable amount of numerical 
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Fic. 2. Variation of wave number with normalized plate thickness for 
complementary waves. Antisymmetric case, pp = V3. 
computation, but if the signal is such that its Fourier representation is 
confined to a narrow band of the frequency spectrum, the phase velocity 
may be represented in terms of frequency over this range by the first two 
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Fic. 4. Variation of wave number with normalized plate thickness for 
mplementary waves. Antisymmetric case, pu 2. 
is terms of a Taylor series, and the group velocity is then given by the formula 
ity | V 


l —y, 33 
wi 1—(w/V)(dV /dw) sis 
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where U’ is the group velocity, V is the phase velocity, and w is the angular 
frequency. This formula cannot be applied in the neighbourhood of a 
cut-off frequency, where the Taylor-series representation ceases to be valid. 
Insuch a region the pulse shape would be greatly distorted, and the concept 
of group velocity would no longer be meaningful. 


Now if a, is the thickness of the plate in unnormalized units, we have 


a ky a; 
ind hence w = kV, = aV,/ay. 
Using equation (47) we obtain 
Ww ap dV Ao dp 
V Ay dw p? da’ 
ind l — Ve (59) 


l+(a/p)(dp/da)  d(ap)/da’ 
[he equations relating p and a are of the form 
d(p) = u&,(p,a), (60) 


: ip? /(p?—p2),/(p2—1 
where d(p) Byvui B (Pp | J (61) 


(2p* pe”)? 


SH(p)Ci(p) 


us, (p, a) 7 for symmetric waves, (62) 
Sa(Pp)Ca(p) 
S1(9)C#(p) : ‘ : 
nd b(p, a) a\Pi alt for antisymmetric waves. (63) 
- p)é ‘(p) . 
in differentiating we obtain 
d ) ous 2 
} n/é é (64) 


da dd dp ous, op 
We deduce the following formulae from equations (61), (62), and (63): 


db 8(1—p?)p?+4y?(3u?—1)p?—8u4p 


, , ~ = - (65) 
ap (2p° pL )°a/( pm fm) (p* 1) 
Cus. )2 2 (1 2 Su 
[y(P°— BCP) VPP — SEP) (66) 
ca Sli p)C#(p)| CH(p) S!(p) | 
Cys, ap | C1(p) : S#(p) | (67) 
pS p)C(p)\Vp>—n2)CK(p) _ (p®— SUP)" 
Cub, | | (./(p?—p? 2)8} (p) \ (p* NCUP)\. (68) 
Ca ( a P)SE(P)| S¥(p) C1\(p) | 
Cub, ap } =P Ch(p) | (69) 
Op Ch(p)S#(p)\./(p?—p?)S#(p) \(p? IC 1(p) yt 
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Values of the group velocity U may be calculated once equations (48) and 
(49) have been solved for p. 

A curve showing the variation of group velocity with plate thickness for 
the principal symmetric mode in a plate having a Poisson ratio of (29 
(corresponding to pu 1-839) is given by H. Kolsky (10), and a set of curves 
showing the phase and group velocities of the first five symmetric and 


antisymmetric modes in aluminium sheet as functions of the product of 


frequency and plate thickness is given by Firestone (11). In the next 
section of the present paper a more detailed analysis is given of a single 
case, and both phase and group velocities are evaluated for all propagating 
modes, as well as the amplitudes of the waves at a point distant from 
the source. 


6. Investigation of the modes generated in a plate whose thickness 

is of the order of a wavelength 

To illustrate the way in which the theory developed in the preceding 
sections may be applied to a practical case, we consider the normal surface 
displacement due to a small circular source of radius 7, on one side of the 
plate. We shall assume the stress on the portion of the surface under the 
source to be unity, and the plate to have a thickness of 2A,/7 and a 
Poisson’s ratio of 4, so that a = k,A,/7 = 2 and p = 2. 

The boundary conditions on the free surface are therefore: 


2=a: 2=0 (forallr), 22 l(@<re), 22-0 (r >F,), (70 


— _ 


Zz a: 2F = 22 = 0 (for all r). (71) 


If we denote the symmetric and antisymmetric components of stress by 
220 and 22°) respectively, it follows that 
231) — 33(2) 4 (r < 9), (72 


5x1) Sa(2 


B™%=—0 (r > Ff), (73 


and from equations (12), (29), and (33) we find for the transforms of the 


stress components y 
ro A (ro) 


2¢ 
g(t) = gall) = 0. (73 


fs() = fal) 


We therefore derive from equations (40)—(43) the following expressions for 


the normal displacement components on the surface of the plate: 


Symmetric component: 





au) — - wna [ J(¢2—1)S# gi Aa(Sro) Ho (Er) a (76 
7 44, 


A a a D,(€;. a, 2) 





; (74) 
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Antisymmetric component: 
a 
44 « 


For small values of 7, 


then become 








IN PLATES 


Ji(r9) Hy? (f 


r) de 


u ro | “¢2@—1)C# C12 (77) 


DAC: a, y) 


we may put J,(C7)) ~ 3¢75, and equations (76) and 


on ieee ae me 
UD) mw alls V(e? 1) su gi to (Sr)é qt (78) 
8Cq4 D(C: a, p) 
2.2 Qyiyr\r dr 
2) am ® f ge—ayog cre Os (79) 
an ° D,(C: a, w) 
don evaluating the integrals we obtain 
Par NT hed “ (1)(/ pit, 
~ —TE ed Vp?—1Se(p)S (p) 28 a LB (80) 
ic,, <— D,(p: 4a, 1) 
1 “ 4 wD itty 
7 TUTTO XN (p2—1)C#(p)C1(y mG PIP ($1) 
1e,, om pP ) a\h a Pp) Di (p: a. pu)” 


the primes denoting differentiation with re 


‘spect to p. 


For large values of 7 


e Hankel function may be represented by the first term of its asymptotic 


pansion, namely o\2 
H(z) ~ (=)*et 
we then find for the displacement components 
n/( 277) r2 ei 4 , * 
~~ Vi Pl p- 
4¢ 14% yi” 
(2a) 2 ye m 
f ap : ~ 5 \{ p(p* | 


$C 44N) a 
we introduce the notation 


D,(p ; 


D,(p) 2ap(p*— 1)}- 
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become 


82) and (83 





2a p(p* 1)! 
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ipr 
1)' Se SIG 5 (S82) 
\P ) oe: a, 2) 
: ip? 
)}CH(p)CUp) . (83) 
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pL ‘1 
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Inspection of Figs. 3 and 4 shows that three symmetric and two anti. 
symmetric waves are propagated, and corresponding values of the phase 
velocity, group velocity, and relative amplitudes of the waves are given 
in Table 2. 








PABLE 2 
| | dpida | viv, | Uv, | Delps22 ®, ( 
n | m Pnm dp/da | = , n(P 3 252) | n(P) 
I I 2°0665 | 01645 | 0°4839 o4175 | 8576 0°01748 
I 2 1°0361 o1659 | O-9651 o*7310 130°3 0°00082071 
I 3 "5050 0°7331 1*9502 05073 | 14°251 O°O71012 
2 I 2°1549 0°1739 0°4577 0*3949 13590 0°03856 
2 2 1°2765 2°7211 0°7532 01488 102°0 o-05610 








The ratio of phase velocity to group velocity for a given wave is a measure 
of the degree of dispersion of the wave, and it will be seen from the table 
above that the values of this ratio for the (1,3) and (2, 2) waves are approxi- 
mately four and five. This well illustrates the need for calculating the 
group velocity in determining the travel time of a pulse of energy along the 
plate and, conversely, provides a warning against facile interpretation of 
travel time measurements in elasticity determination. 
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SUMMARY 


Second-order perturbation theory, using known solutions for a circular cylinder, 

3 applied to the exact elasticity equations to determine the wavelength frequency 

ives for the longitudinal vibrations of a non-circular cylinder. General formulae 
given and applied to elliptic, square, and rectangular beams. Owing to 

sonance, the simplest technique is inadequate at certain wavelengths for elliptical 

nd rectangular beams, and this makes the method rather cumbersome. It is also 

ear that such sections as a square deviate sufficiently from a circle as to demand 
inclusion of third-order terms 


1. Introduction 
fue velocity of longitudinal stress waves along uniform beams has been 
culated for sections which are circular and rectangular using the exact 

sticity equations, assuming Hooke’s law for an isotropic material (cf. 
R. M. Davies (1), for references). In this paper and a later one we seek 
the solution of the same problem for a beam with a general section whose 
boundary, in polar coordinates, has the form 

a(l > b.cos s@), (1.1) 

vhere the coefficients b, are small. We do this using a perturbation theory 
used on the solutions in polar coordinates used to solve the problem for 
. circular cylinder r = a. In other words, for a given wavelength, we 
expand the phase velocity as a power series in the coefficients b,. 

In this way we hoped not only to obtain more general results than those 
btained previously and to observe the effect on the dispersion curve of 
ious changes in the shape of the cross-section, but also to explore the 
possibilities of perturbation theory in this type of problem and to see if it 
ld be used to derive other approximations which might give better 
esults. Moreovet by using a series of type (1.1) as an approximation to 
sections with sharp corners, it may be possible to estimate the effect on the 
elo ity of rubbing off the corners. 

The construction of a second-order perturbation theory was straight- 
lorward, but an extension to third-order is very elaborate. It is easy to 
nclude what seem to be the main third-order terms, but the magnitude of 
the residue could not be determined. This meant that the limits of validity 
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of the second-order theory, which fixes bounds on the magnitude of the 
coefficients b,, could not be determined by an examination of higher-order 
terms. We therefore apply our results to a square cross-section, for which 
experimental results are available (2). Although this comparison was 
rather inadequate, it made two things clear. First, that the second-order 
theory is not adequate for such a section, and higher-order terms should ly 
included, a result not unexpected, as a square is appreciably different froma 
circle. Secondly, the effect of the higher coefficients b, with s > 8 was quite 
negligible, i.e. less than 1 in 10°, for all frequencies corresponding to velocities 
greater than the shear wave velocity, and possibly for higher frequencies, 

As the work is in part exploratory, and the calculations lengthy, they 
were not made for many frequencies, and the frequencies chosen wer 
near crucial points in the dispersion curve. For example, a value giving a 
phase velocity approximately v2 times the shear wave velocity was chosen 
because of a limitation on the perturbation theory which was not suspected 
when the calculations began. This limitation is due to a degeneracy of two 
distinct modes of vibration of a circular cylinder in this region, and needs 
a little explanation. 

As the basis of our theory we use solutions of the elasticity equations in 
polar coordinates, which can be classified according to the number (n) of 
nodal planes intersecting the axis. The fundamental modes of vibration 
of the circular cylinder are then classified by the solutions needed to describe 
them. For the longitudinal and torsional modes n = 0, for the flexural 
modes n 1, for the screw modes n 2, and so on (the existence of higher 
modes is not yet certain). Using an approximate theory it was discovered 
that the longitudinal and another mode, which was later described as a 
screw mode, had dispersion curves which intersect in the region c? = 2c’, 
and that the torsional and screw modes intersect where the latter crosses 
c? = c?. Now the perturbation theory breaks down whenever the boundary 
conditions couple two solutions which describe degenerate modes of the 
circular cylinder. Hence it breaks down for the present calculation of the 
longitudinal mode near c? = 2c? and gives poor results in the neighbour- 
hood of this value. 

This difficulty does not arisc for a square cylinder, where the solutions 
needed are those for n = 0, 4, 8,... and there is no coupling between the 
longitudinal and screw modes, nor does it arise for sections such as al 
equilateral triangle. It does arise for an elliptical or rectangular cross- 
section, where values of n differing by 2 are needed, and as a result, the 
dispersion curves for the longitudinal and screw modes, n = 0 and 2 


respectively, are replaced by two adjacent branches which do not intersect. 
This is in agreement with the experimental results obtained by Morse. 
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The results of our calculations are summarized in section 3 and Table 1. 
Using this table it is possible to calculate the change in phase velocity at 
various wavelengths, due to changes in the cross-section due to terms 6,, 


8 6). 


2. Method of calculation 

To find the oscillations of a circular cylinder it is appropriate to use 
cylindrical coordinates (r, ,z) and to set up the elasticity equations for the 
components (w,v,w) of the displacement in the radial, transverse, and 
ongitudinal directions, respectively. 

The axis of the cylinder is along the z-axis, and we look for solutions 
where the components are proportional to e4”9+*:-«), where n is an integer, 

27 A) the wave number, w (= 2zv) the angular frequency, and ¢ (= w/k) 

the phase velocity. Thus ka is the number of wavelengths in the perimeter 
of the cylinde1 

There are three solutions of the equations for each value of n, which are 


labelled 1 l. z= 3 


u u,,(A,,;, cos né+ B, , sin n8)cos(kz—wt) 


v = v,,(A,,;sinné— B,, cos n@)cos(kz—wt) }. (2.1) 


Ww w,(A,,; cos n6+- B, , sin n@)sin(kz— wt) 


ni n 


The expressions for w,,;, U,;, W,; are given in the Appendix. The general 
solution is a linear combination of these solutions, and the stress compo- 
nents Dow, etc., are determined by the formulae given by Love (3). 


In the case of a circular cylinder, the boundary conditions, 


Pz, Pro Di. 0 atr=a, 


j 


we satisfied by combinations of the solutions for each value of n. In fact, 
we only need the solution n = 0, ¢ = 1 for the torsional mode and its 
harmonies, and a mixture of the solutions n = 0,7 = 2 and 3, for the longi- 
tudinal mode of vibration and its harmonics. The two flexural oscillations 
A, and B, are determined by the solutions with n = 1, one needing the 
solutions with coefficients A (ef. (2.1)) and the other solutions with the 
coefficients B. The two screw vibrations A, and B, require solutions with 
2, and so on. The equations which determine the phase velocities are 
the Pochhammer equations, which have been solved numerically by R. M. 
Davies (4), and Bancroft (5) for n = 0, and by Hudson (6) for n = 1. 

By plotting phase velocity against wavelength or frequency for each 
mode, a branch of the dispersion curve is obtained. These branches and 
their relative positions are shown in Fig. 1, for long waves. This diagram, 
which does not include the harmonics, shows that the number of inter- 

. 
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sections, or cross-overs, of the branches is small in this range of wavelengths, 
There may be more when the wavelength is small compared with the mean 
diameter. These cross-overs are important because any small change in 


the boundary which mixes two modes of vibration, causes a separation of } 


their dispersion curves at the cross-over, as indicated by the dotted lines, at 














some points in the figure, which cannot be dealt with by simple perturbation 
theory. In fact, our calculations break down near such points, although 
this type of degeneracy can be dealt with easily in principle, because the 
amount of numerical work is excessive. If we include harmonies of the 
different modes there are further cross-overs with the screw branch and its 
harmonics. 

As we see later, the solutions for a normal oscillation with the perturbed 
boundary (1.1), which only involves terms in cosines and therefore has one 
plane of symmetry, divide into two groups containing, respectively, the 
solutions (2.1) with coefficients A and with coefficients B. The first 
includes the longitudinal mode and the second the torsional mode. Detailed 
considerations depend on the symmetry of the cross-section and the 
coefficients b, which occur in its equation. In general a non-zero coefficient 
b, mixes the elementary solutions with values of n differing by a multiple ofs. 

For example, a square boundary only requires terms in ,, bg,..... Conse- 
(0, 4, 8,...), (1, 5, 9, 
(2,6,...), and (3,7,...). Thus, not only are the torsional and longitudinal 


quently the solutions are mixed in the four groups n - 


vibrations, n 0, never mixed, as explained above, but they never mix 


Thus 


1 or the screw vibrations 2. 


with the flexural vibrations n 
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hs the cross-overs with the screw vibrations for the circular cylinder are 
"aN | unimportant for our calculations, but are retained in the dispersion curves 
“| foranexactly square beam. The same is true for a beam whose section is 
101 ? an equilateral triangle, whose motion is described by mixtures of solutions 


40 | with values of n differing by multiples of three. 

This is not true for a rectangular or elliptical beam, where terms in b, occur. 
There are now only two groups of solution, n even or n odd, for either A 
r B,so that the longitudinal and screw mode A are mixed, and the torsional 
ind screw mode B (7). The flexural mode (n 1) is mixed with higher 
rders. The low-frequency cross-overs are now relevant, since the separation 
of the dispersion curves means that the low-frequency part of the longitu- 
dinal branch joins up smoothly with the high-frequency screw branch, and 
vice versa, together with a similar effect on the torsion branch. 


The calculations for non-circular beams proceed as follows. The boundary 





onditions for a uniform infinite beam, on the curved surface, are 

p,,r de pjgdr Q), (2.2) 
where the suffix j) has the values 1, 2, 3 corresponding to suffixes r, 0, z. 
Both r and the stress components on the boundary are periodic functions 
ff @ with period 27, so that these equations can be expanded as Fourier 


| series, and become 


> (F;,,, cos mé+- F;,sinmé) = 0 (9 1, 2, 3) 
F101 m 
ugh nf 
the where 7K. | cos mO{rp;.—(dr dO) p jos dd (m #0) 
the | <f (2.3) 
1 its nF ( (rp (d) dO) p 9} dé (m 0) 
0 


bed | The equations for F%,, are similar to those for F,,,. Then the orthogonal 


jm jm 


one property of the circular functions enables us to replace these by the equations 

the 

first Fin U I jm V (all m, )). (2.4) 
.4 | Now the stress components and the Fourier coefficients are linear functions 


if the displacements and the coefficients A;,, B 


n? in 


in the displacements. 


——EE 


7 Moreover, it is easily shown that the coefficients B do not appear in F,,, 
of s nor the coefficients A in F’.,, when r is an even function of # on the boundary. 
ii Hence these equations form two sets of linear simultaneous equations 
iY din pv’ . j , . of 
F > Pir A,, 0, Fim = > Qj", B;, = 0 (allj,m) (2.5) 
lina ‘ = 


mix | Which have non-zero solutions only when the determinants of the coefficients 


‘hus vanish, 


Pe 0, qin 0. (2.6) 


jm 
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The determinant P includes the solutions n = 0, i = 2 and 3 which 
appear in the longitudinal vibration of a circular cylinder, so that it has 
one solution corresponding to the longitudinal oscillation of the given beam, 
The determinant Q includes the solution n = 0, i = 1 for the torsional] 
oscillation of a circular cylinder, so that it has one solution corresponding 
to a torsional oscillation of the given beam. We can also deduce from 
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Fic. 2 
equations (2.3) and (2.5) the earlier statement on the mixing of the various 
solutions, for example that, when y contains only terms in 26, 44,..., then 
the elements P}"., Qi" of the determinants are zero unless m and n are both 
even or both odd, so that the determinants factorize into two smaller 
determinants. 

The solution of (2.6) by perturbation theory begins with the expansion 
of the elements in powers of the small quantities b,, which we may assume 
to be of the same order of magnitude. The determinant P has the form 
shown in Fig. 2. 

When the &, are all zero, this determinant reduces to a product of small 
determinants D, D, D,..., where D, is constructed from elements with 
m = n = Oand has two rows and columns, D, is constructed from elements 
m = n(n > 0) and has three rows and columns. When the b, are not zero, 
there are further terms of order b, or higher order in the off-diagonal 
positions D,,,, (m ~ n). 
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For a complete solution, the determinant should first be transformed so 
that all the off-diagonal elements D 


ban (m ~ n) become zero. In particular, 
to find the longitudinal vibration, we should transform it so that the terms 
D,. and D 


0 m{ 


one corresponding to ), which gives the longitudinal frequencies, and the 


, are zero, whereupon the determinant factorizes into two parts, 


ther including all remaining motions. Such a solution is not possible, but 
tis possible to obtain answers correct to the second order by removing the 


terms in 6b, from the elements D, 


bn: This can be done fairly easily, by 


subtracting from these two rows, multiples of the three rows of D,, i.e. by 


subtracting multiples p, q, 7 of the three rows m = n, chosen so that 
rin din iN » arin , » 2 
in Pn I 2n 4 Tan % a 10 (2 l, =, 9) ) (2.7) 
: ; ° of 
P= Pn Pon q P3n rn 20 (2 ? 2. 3) 


In these equations the terms P‘” 
| 


‘* are needed only to order zero, so that 
they have solutions, unless D 0, without any limitations on the right- 

ind sides. The subtraction is performed for every value of n (except zero), 
are of order 


with the result that the terms in the position of the matrix J), 


ind the determinant LD, has become, with revised elements, 
pio pio 
’ 10 20 9 
Do _|* (2.8) 
220 P21 
PX PS 
Since the elements D,,,, (mm 0) are at most of order b,, the whole determi- 


nant can be expanded in the form 


Dy X {D, Piss terms due to b,!+ terms of order 63 0. 


j 


Provided the determinants D,, D,,... are not too small, this equation has 


me solutior 
alii ” dD 0 (to order b?). 


0 


In the evaluation of the correction terms we assume that the phase 
velocity differs little from the solution of D, = 0 for the longitudinal 
mode of a circular cylinder of the same wavelength. Our perturbation 
theory therefore breaks down either for large values of the coefficients ),, 
t when one of the determinants D, becomes small for this velocity and 
wavelength, i.e. the two dispersion curves come close together, for two 
lifferent branches of a circular cylinder corresponding to zero values of 
Dand D,. In fact, for wavelengths which are not much smaller than the 
ateral dimensions of the beam, only one determinant becomes small in 
this way, namely D,. This does not matter in calculations for beams with 
lourfold symmetry such as a square, as explained above, but it does upset 
those for a rectangle or ellipse where n = 2 is also required. For such 


beams calculations by our present method are inaccurate in the range of 
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wavelengths where D, is zero or small. Points are obtained for both smaller 
and larger wavelengths, but the former are probably not very good, and 
although they describe essentially longitudinal vibrations, they do not 
lie on the same (continuous) branch of the dispersion curve as the latter 

The various elements of the determinant used in these calculations are 
given in the Appendix. 


3. Results 

In view of the numerical work involved, the calculations were made fo 
only one value of Poisson’s ratio, co = 0-3, for the limit of long wavelengths 
and four other frequencies, one of them the frequency such that the phase 
velocity c equals the shear wave velocity. The results are summarized in 
Table 1. 


The final determinant Dj is expanded and expressed in the form 
—Dy = A+ 


> 362.0,, (3.1 


where 


A = {1+¢9(Ba)} 
K 1+A(k?+- «?) 2a", 


4(1—B?/k®)(1+-K,(2a)). 
box) 


The coefficients [, are functions of the wavelength. The phase velocity is 


and 


xJy(x)/Jo(2). 


now determined from the equation Dj = 0, or 
A > 3462 T,. (3.2 

The phase velocity can be expressed as a function of A in the form 

c7{ AA+ BA?+ CA?-+...}. (3.3) 


c? (section) —c? (circle) 


The coefficients A, B, C are givenin Table |. This table includes, therefore, 


TABLE | 


Values of VT, 

















(a) (b) (c) (d) (e) 
c*/cX(circle) | 2°6 49/19 7/3 2 i 
Ba ro) 052614 1°53263 184118 ° 
xa/1 ° 0°21450 0°76031 1°20534 3°6232 
I, ° 0°82874 28-3785 12*7090 45°529 
Re ° 0°78093 10°9652 19°2997 4°430 
I, ° 0°78416 11°3056 226800 | 10281 
i, ° 0°79580 12°5337 26°6755 | 19°383 
Tio ° 140881 31°379 30°225 
Lis °o os 1O*1O5 36° 386 42°109 
A ° 0°29299 0°51133 O°91139 0°5 303 
Bb | © | 10000 0*3900 0°86338 O°21245 
O15 


Cc a: ie ia ae fe) 
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uller | all the information necessary to obtain a number of points on the longi- 
an tudinal branches of the dispersion curve. 

not Using these results, phase velocities were calculated for three sections, 
ster in ellipse of eccentricity « }, a square, and a rectangle with sides in the 
are | ratio6:5. These were chosen for reasons which appear later. The coefficients 


b, for these sections are given in the Appendix, together with the radius a 
of the cylinder with which they are compared. This radius a is the mean 


value of the radius vector r around the section. 











1 th 
PABLE 2 
rt 
L/P; ingular frequency; a mean radius. 
last 
wa Ellipse Square | Rectangle 
6 - - 
2°5786 | 2°5777 ‘577! 
) 2°0275 3 2°3140 2°3016 2°2549 
2°01 46 1°8921 1*9450 
$°299 I 10304 0°9834 10558 
The results are given in Table 2. In judging the results it should be 
Vis remembered that the coefficients 6, in all three cases are probably as large 
is can be justified by using second-order perturbation theory. Moreover, 
) (3.1) shows that, when c? = 2c and 8B = k, the value of Ba is completely 
letermined for a circular cylinder, and is independent of Poisson’s ratio. 
This value of the phase velocity, according to the approximate theory, also 
makes D, vanish as well as D,, and although the true value for this to occur 
is a little less, it is close enough to upset the perturbation theory, for a 
re . 


rectangular or elliptical section. In fact, we see from Table 2 that for the 
? ellipse the points for lower frequencies lie below those for the circle, and 
lie on one branch, and those for higher frequencies lie above those for the 
circle and lie on another branch. A similar comparison between the square 
ind rectangle gives the same result, for sections with the same value of a. 

We can compare these results for a square and rectangle with those of 
Morse (2) assuming, as he does, that the shear wave velocity is 2,160 metres 
per second and Poisson’s ratio is 0-30. This enables us to calculate the 
lrequencies corresponding to the three phase velocities near c?/c? yf a 
ind 1-5 (the set for 1-5 is estimated from the calculated values). The 
results are given in Table 3.+ The difference in the values for a square at 
the lowest frequency can be explained by an error in Poisson’s ratio, but 


this does not explain the difference at c? = 2c?. This can only be due to 


We are grateful to Dr. Morse for his actual values which enabled this comparison to 
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an error in the shear wave velocity of 3 per cent., which we can rule out, 
or to an error due to the use of perturbation theory. The error in the values 
for a rectangle may well be due to the fact that we are close to a zero of D,, 
The fact that the difference between the calculated values for a square and 
a rectangle is less than the experimental may have the same cause. 

In this table we have included for comparison calculated frequencies for 
a circular rod with the same cross-section as that of the square, both for 


G 0-30 and 0-35, the latter value being chosen by Morse for this purpose 























mh ine 
TABLE 3 
“/ 7) J or , ( le 
( Square Rectangle ircle 
= ———__—— ~ - —- 
C3 (exp.) | (calc.) (exp.) (calc.) \o o30\a 0°35 
(c) 2°03 35350 3,270 3,260 3,250 35295 39317 
(d) | 2°60 3,050 2,980 3,240 3,040 3,053 3,053 
Est 3700 2,700 2,680 2,810 2,740 2,750 2,751 








In conclusion, the perturbation theory developed does not seem to be 
accurate enough for a square, using the circle as the unperturbed shape. 
but more experimental results: for various sections are really necessary. 
One positive result, however, does emerge. It is that the two branches 
observed by Morse agree very well with the two branches obtained by an 
interaction between the longitudinal and screw vibrations, and are not the 
longitudinal branch and a harmonic. 
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APPENDIX 


Apart from the factors involving z, ¢, and 6, the solutions used in section 2 are: 


d ] D 2 D 3 
Uni nJ,(Br)/Br J, (Br) Jy(ar) 
Uni J) (Br) nJ,(Br) Br nd, (ar) xr 
W yj 0 B.J,(Br)/k kJ, (ar) /x 
where k = 277/A is the wave-number, and « and B are defined by 


pw*/p B? +k, po? /(A 2) a? -+ k, 
The dash denotes a derivative with respect to the argument Br or ar. These include 
n 0. When n 0, u w 0 for 7 l, and v 0 for? 2, 3. For the solutions 


i 1, 2 the dilatation 6 0 and for i 3,6 (x? +-k?)J,, (ar) /a. 
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ut The following are the elements of the determinant P, to the desired order in the 
” efficients 6,, using 2 Ba and y va. (K is defined below equation (3.1).) 

D n. P Ji(a) + DY 4b2{x?2Jy(x) — Jo(x) — 26 (x)/2} 
P (B2— k?)/2Bk[ Ji (x) > pbrz{axJo(x)-+4 x? J} (x)}] 
‘ p® — J(y)—(K—1)A(y) + ¥ POA(y?— Joly) — 2S o(y)/y —(K — 1 (yJo—y*p)} 
P kJi(y)/a + > 4b kat Jo(y) + yJo(y)}- 
10 D.: P i ] i 2 i 3 
f n(aJ’ —J,)/2 4 J” (x) MIi(y)—(K—1)J,(y)} 
” 2? J ha? J n*J,,)/ 2a bn(J—aJ> )/x? n(J,, yJ,) 2y? 
} nkaJ, /42 (k?—B?)J* /4kB kJ) (y)/2a 
- P¥ 
nb, {J,(1—}a rd) \ |x 0 
2 4b, {xJ)(1 n? J, \/2? bb, ka(k? — B?)J,,(ax)/k* 
3 lb j n?)J,,/y2+ J; /y—K(J,+ yJ,,)} bb, kaJ,(y) 
P Pye 
b, Jo(x)(1—2?)/a bb, {Jo(y) + J6/y—K(Jg+ yJ0)} 
2 snb,, J lnb, {Jo + Jo /y—KA(y)} 
3 Lb (k? ; uJ, r)/k? 4b, kaJ,(y) 
} ficients rious sections. The following values were used: 
l lipse of eccent ty ¢ 1. 
ee b, 0-07184, b, 0-00387, b, 0-00023, b. 0-00002. 
Square of side 2s, where a 1-128: 
b, 0-13941, b, 0-04397. 
Rectangle, witl des in ratio 6:5: a -O218d, where d is the longe r side: 
0-10985, b, 0-12713, b, 0-04359, b, 0-03148. 
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SUMMARY 

The second order perturbation method developed in the previous paper (1) js 
adopted for the consideration of torsional waves in a bar of general section. Thi 
determinantal frequency equation is obtained correct to second order, for the first 
mode of torsional oscillations of the bar. This equation is solved for particular 
wavelengths and dispersion curves sketched for bars whose sections approximatt 
to a square, rectangle, and ellipse. As in the longitudinal case it is found that th. 
method breaks down for the rectangle and the ellipse at a certain point due to 


degeneracy between the torsional and screw modes. 


1. Introduction 
[nN the previous paper [ the theory of longitudinal vibrations of an infinitely 
long bar of elastic material having the cross-section 

r = a(1 + > b, cos 86) (1.1) 
has been examined by means of a perturbation method. The method is 
equally applicable to the case of torsional vibrations which will be con- 
sidered here. 

The torsional vibrations of a bar of circular cross-section have been 
investigated by Chree (2), Bancroft (3), and Davies and Owen (4). In 
this paper the results obtained for the bar of general section (1.1) are 
applied to cylinders of square, elliptical, and rectangular cross-sections 
whose Fourier expansions are terminated at convenient points, to obtain 
the appropriate dispersion curves for the fundamental mode. The results 
obtained being correct to second order in the coefficients b,. 

As in the longitudinal case it is found that at certain wavelengths 
degeneracy occurs between the fundamental mode of torsional vibrations 
and the screw mode. This does not affect the perturbation theory for a 
square, as pointed out in I, but the method is not valid in this region fot 
the ellipse and rectangle. 

2. Method 

The method employed is exactly similar to that of paper I. The general 
solutions of the elastic equations are obtained in polar coordinates and 
made to satisfy the boundary conditions for the bar. On expanding the 


[Quart. Journ. Mech, and Applied Math., Vol. X, Pt. 1 (1957)] 








stresses 


simulta 
of the 
solutio 


These | 
P' col 
the tor 

For 
where 
here | ¢ 
every 

The 
equati 
by Do 
all ele 
As in t 
b?, anc 
from t 
first 1 
for th 
bs, by 


where 

We 
is asst 
the te 


spond 


3. Re 

The 
polar 
eleme 


Th 


wher 
and 4 


Fo 

















—_ ws 











VIBRATIONS OF BEAMS. Il. TORSIONAL MODES 


‘ 


stresses in these boundary conditions in Fourier series a set of linear 
simultaneous equations in the arbitrary constants of the general solutions 
f the displacements are obtained and the condition that these have a 


solution gives rise to two determinantal equations 


i,j 1, 2, 3, 
PRi=@; i@ei=e ( - 
m,n DB Bios 


[hese determinants have been discussed in paper I where it is shown that 


(2.1) 


P corresponds to the longitudinal oscillations of the cylinder and |Q) to 
the torsional oscillations. 


For a circular cylinder @ consists only of diagonal minors D,, 


BA, Bhyn, 
where J, is (1x 1) and PD, is (8x3). For the general cylinder considered 
here ( contains the same terms together with terms of order b, or more 
everywhere 
The modes of torsional vibrations of a circular cylinder are given by the 
equation D, — 0. Similarly, those for the general cylinder will be given 
by D, = 0, where Dj is the new first element of |Q| obtained by reducing 
elements in the first row and column, apart from the first one, to zero. 
\s in the longitudinal case, however, we only consider the solution to order 
ind it is therefore sufficient to eliminate from the first row of |Q!, apart 
from the first term, all terms of order b,. This is done by adding to the 
first row suitable multiples of all the other rows. The frequency equation 
for the torsional oscillations of the general cylinder is then given, to order 
by the equation D’ 0. (2.2) 
where Di is the resulting element in the first row and first column of |Q). 
We will consider the first mode of vibration only. As the perturbation 
s assumed to be small, the multiples of the rows of D, which eliminate 
‘terms of order b,, in ),, are evaluated at any point by using the corre- 


sponding solutions for the circular cylinder at that point. 


3. Results 

rhe general displacements obtained by solving the elastic equations in 
polar coordinates are given in paper I (Appendix). These give rise to the 
elements Y:" shown in Table 1. 

lhe frequency equation for torsional vibrations of a circular cylinder is 

D, 1 +34 (x) 0, (3.1) 

where d, (x) xJ,(x) J (x) 

nd x is defined in paper I. 

For the fundamental mode of vibration this has the solution 


0, 
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Using this value of x to eliminate the terms of order 6, in the first row of 


(| we obtain the frequency equatic 
order L?, : 
Ds 


where (1/I.) 


1 + (2s—3)z?/2s(s? 


and z and @, are defined in Table 1. 


i+ sho r 


m. for the general cylinder, correct to 
> 422sb?T, = 0 


1)—8. 


s 


, (3.2) 


24 {45?(s +] (1 -6,)} 


TABLE | 


Elements Q'”, of 


the determinant Q 











. 10 € Y 293 9 2 6 
Dy: Qin = 14 45/2+ > 03 {1+ o/2—x?(1+ do)/4}/2 
D Qn i l i 2 0 3 
) ] 1— 26, /n+n@,, n—6, ] Ky6, n+né 
) 2 9 5 
} «= n 9, l ween aN nO, m On 
) 3 (1—a?/z?)/2 0,,/2 l 
l 2 
Don: Qi = 5,(8,—278, —n) 
8 Fy b,(1 2x? 4 a 26 »n nO.) 
Qi5 b,(@, — Ky*6, —n) 
10 9) /6 
D 0 im Dn m(1 T do 2) 2 
Lo 9 2/9 
Qom Bin ( 1 — }o/2+2x?/2) 
10 
Q3m b ni—m 4) 
where do xd o(x)/ Jo (x) 
9, nJ ,(x)/xJ (x) 
6, nJ,(y)/y J? (y) 
2 ka 
ry » 
TABLE 2 
y aw . > . 
Values of i. (s 2, 4,6 8). o 0-3 
(ka)? Al “ re + & 
inlmenataan — . |- _ 
° 1*0000000 I*0000000 I°0000000 | 1*0000000 
I*4 1'0526623 o'9gg9T1516 0°9934799 | 9°9956126 
2°8 I°1193416 0°98 32335 0°9867866 0°9913374 
$°2 1*2041572 | 0°9761920 | 0°9812849 | 0°9871714 
5°6 1°3133640 | 0°9699798 | 0°9755883 | o-9831112 
77O +) 1°4568390 | 0°9645553 | 0°9701457 | 0°9791562 
}° 3755101 


The values of [, for s = 2, 4, 6, 8, 


values of ka and a Poisson’s ratio o 


0°9451 332 


0°9464404 | o-9608491 


are given in Table 2 for a number of 
0-3. This enables the phase velocity 


at these wavelengths to be determined for any bar whose cross-section can 
be put in the form (1.1) provided the }, are sufficiently small. Equation 


(3.2) has been solved to give these phase velocities for bars of square, 


rectangular, 
taken up to s 


and elliptical cross-sections whose Fourier expansions are 


8, and the resulting dispersion curves are shown in Fig. |. 
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VIBRATIONS OF BEAMS. Il. TORSIONAL MODES 


4, Degeneracy 

\s has been pointed out in [, the method applied here breaks down 
when D, and any PD, vanish together. That occurs when the dispersion 
urves for the torsional mode of vibrations of the circular cylinder and 
any other mode intersect. It is in fact found to happen for ka ~ 4, when 
the torsion curve intersects that of the screw mode D, = 0. However, 

















L ; i 
1 2 3 ka 4 
Fic. 1. Dispersion curves for torsional vibrations 
a) Cireular cylinder (b) Ellipse (eccentricity }) 
Square (d) Rectangle (side ratio 6:5) 


is explained in I, this does not affect the results for sections such as a 
square or triangle which do not contain a term in b,, but will affect the 
results for a rectangle and ellipse as can be seen from Fig. 1. 

Calculations have been started to obtain the dispersion curves for the 
screw vibration D, = 0, for different values of Poisson’s ratio. A sketch 


curve for o 0-3 through the three points given in Table 3 is shown in 


Fig, 2 


TABLE 3 


Points on disp SiON Curve for screw vibrations (o = 0-3) 
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4-0r 











% | 2 3 ka 4 





Fic. 2. Dispersion curves for (a) serew and (b) torsional vibrations 


of a circular cylinder 
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Received 23 February 1956] 


SUMMARY 


[he problem of a magnetic dipole acting as a transient current source when 
tuated over a two-layer earth is considered by using the methods of integral 


transforms. The general expressions have been applied to deduce the induced field 


he surface of the earth when the dipole is located on the surface, in the cases 
rhen the earth has homogeneous conductivity, and when the conductivities of the 


two layers are nearly equal. The expressions obtained for the induced magnetic 


1 components are in a form suitable for numerical calculation. 


1, Introduction 

Usrna some results obtained by Price (1), Gordon (2) has considered the 

problem of an oscillating magnetic dipole outside a semi-infinite conductor 
nd Bhattacharyya (3) the problem of an oscillating magnetic dipole over 
two-layer earth. 

In this paper, using the methods of integral transforms, we consider the 

ise of a magnetic dipole, acting as a transient current source, when it is 
situated above a two-layer earth. In particular we confine our attention 
toa step function and an impulsive current source situated over a homo- 
geneous earth and also over a two-layer earth where the conductivities of 
the layers are nearly equal. 

Expressions for the induced magnetic field components are deduced in 
forms which can be readily calculated numerically. For geophysical 
exploration this may prove useful for the interpretation of field data. 

The use of integral transforms enables the field components to be obtained 
lirectly from Maxwell’s equations, avoiding the use of a Hertzian vector 


ind the vector and scalar potential functions. 


2. Statement of the problem 
A magnetic dipole is situated above a semi-infinite two-layer earth. 
Cylindrical coordinates (p,¢é,2) will be used in which the z-axis is along 


upward vertical. Let the surface of the earth be taken as the plane 
and the source be located at (0,0,a) in a region 1 defined by 


O< sce 


Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 1 (1957)] 








80 J. S. 


LOWNDES 


with characteristics o,, €;, “,, and let the earth (z < 0) be considered 


stratified horizontally into two regions, 


Region 2, 0 > z h, characteristics, cy, €5, [s, 


Region 3, —h >z 0, characteristics, o3, €3, 15- 


The characteristics o, «, 1 denote respectively the electrical conductivity 
dielectric constant, and permeability. We assume axial symmetry an 
denote the current density in region | by J4(p. z,¢). 
(in m.k.s. units) are then 








Oz ot 
é ' oH,. 
(pig) cea (1 
p Cp ot 
OH,, oH;. ( é 
“4 % = |o,+¢;,—)E,4+J4(p. z, t) 
dz op \' ‘a wT eglp 
where i = 1, 2, 3, denotes the components of the separate regions and 
J, = 0 when i = 2, 3. The boundary conditions are 
(a) H,, = A, (d) H,, = Hs, 
(b) Ex4 Fxg, at z 0, (e) Ey : E34, at.z h 
(Cc) py Ay. = be Ag. (f) 2 Hy. = bs As. 
Applying a Laplace transformation in ¢ defined by 
E.g(p.2,8) = | E,g(p,z. the dt 
0 
to equations (1), we find 
‘ Big hs H,, 
Cz 
l C + ) 
(pL) ~P;zs H,, . (- 
p op 
oH,, 2H,. 5 7 
= © = = (o;- €;S)E-4 +d 4(p. 2.8) 
Cz cp 


Transforming equations (2) by a Hankel transform defined by 


[ pJ, (0p) H;.(p, z, 8) dp, 


0 


H AG, Z,8) 





Maxwell’s equation 





of ordet 


Solving 


where 
The 


be give 


Region 


Sine 


Substit 


and 


Region 


Sine 


Reqi 
Leqor 


Sin 








red 
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order zero for H;. and order one for Jg, H,,, E.4, we have 


dE: 
dz lad SH ; 
dé BK; SH ;.. 
AH 
— Ax (a, €; 8)é i4 + J (9, z 8). (3) 
Solving the above for G4 we find 
dé 9 ¢ 4 
re Sig = His F4(8, 2,8), (4) 
ynere - 0 A; S\o . 


rhe solutions of (4) corresponding to each of the three regions can now 


be given 
Region | 
Since 6 44 >Oasz—> yf 
[1,8 ’ . F ~ 
or = $4(O,A,s)e-ms-A dy 4 A(6, s)e-™ =, (9) 
aad b | 
Substituting (5) into (3). we get 
H i ff] é 
M — | J4(0,A,s)e-m™=- dd Ae-™®, (6) 
= . f,8 
l é Z A 2) y ~ 
: | (-).44(0.a, sem = dx — 1 Ae-n (7) 
: 2 a 2 A P bys 
her VU] 4 
Since /, 0 we have 
6», B(6, sje C'(6, s)e” 
° 7 
Uf /2 Be 1 le (e1 
Mod Bos 
G 6 ' ; 
a Be- Cene-, (8) 
7 [Lo 8 fos 
m3 
since f 4 0, 634 > 0, as 2 2, then 


6 34 D(6, s)e%s- 


1): 
Yy 13 Dens-. 
, 3° 


‘ed 


IBY, 3=, 
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The problem now reduces to the following: given a particular form of J, 
and applying the boundary conditions, to evaluate the terms A, B, C, D, 
and hence the field components. 


3. General solution of the problem 
We regard the dipole, located at z = a, as an infinitesimal loop with its 
axis in the z-direction. The dipole can be represented by the form 


__@ b)d(z—a 
Oe Ala. iin 

b—0 2mp 
where 27p is the normalizing factor and 8() the Dirac delta function. The 
current density is then seen to be 


J4(p, 2, t) lim MT (1) ea) (10) 
b 


2 


~0 7p 
where M is the magnetic moment of the current loop and 7'(f) is a function 
of time only. 
Transforming (10) with respect to¢ by means of a Laplace transformation, 
Js(p,2,8) = lim M Tg) 2lP - _- a) 
b—0 7 p~ 


then applying a Hankel transformation in p, we find 


$ 4(9, 2,8) MT(s)8(z @) ys 21098) 
ed bo O 
Proceeding to the limit, x 
$ 4(9, z, 8) M I (s)00(z—a) an) 


9 
aT 
Substituting (11) in (5), (6), (7) we have for the transformed field com- 


ponents in region 1, 


— MO, sT(s)_ 


6G m2-a@ | Ae m=, 
” 477, 
MET (s 0 
‘a 16 T's), Wwso-a Ae-™?, 
; din; fj, 8 
J Z - 
KH, 10 4 a T'(s)e miz—a\ n Ae mz, (12) 
‘ 4a \|\z—a [18 
Henceforth we set .; = 1, and apply the boundary conditions to (8), (9), 
and (12): we find 
Mé6sT(s 
(a) — . , me A —No B+ C 
Arr a ci 
MésT(s 
(b) nae 1 sT'(s) , m4 _} A B4 ¢( : (13) 
47m, 
(c) No Bet2"-- yn, Ce-m:" = yn, De-'" 
(d) Benzht-Ce-nh De-nsh 











Solvin 


and su 


The fi 
and th 
with tl 


Simila 


where 


The al 


compo 


4. Th 
For 


We 


(ii) 


(ili) 


Sub 


Hanke 


where 


Phe 

















\ TRANSIENT 





MAGNETIC DIPOLE SOURCE 


Solving the above for A we find 


A M6s1"(s)| (2+ 1)(M2— 13 (2 M2 M3) e-m@, (14) 
77) (72— 1)(N2— Hg )e- “= — (Ne-+ 91)(N2+ Ns) 





ind substituting for A in (12) we get 


M6@?T(s) 
a 


, m2—a) | wt. (15) 
71) 


[he first term in the above expression corresponds to the inducing field 


ind the second term, -# f., is the induced field. We are here only concerned 


with the induced field 


M*. ee 1 (1, 2, He. (16) 
i 
Similarly w* a he Ngee t®, (17) 
whee - 
Liss, tae) (2+ M1) (N2— Ns)é eticnme. ~My) 2+ Ns) (18) 


; 272 - . 
(H2—1)(M2— 93 )e-“"™ — (N2+ M1 )(2+ Ns) 
[he above are the general expressions for the transformed magnetic field 
omponents 


4. The induced magnetic field over a homogeneous earth 


For a homogeneous earth 7», = 73. In this case (18) becomes 
Yo 7 
I(n,, 2) = 2—". (19) 
N27 1 


We will now make the following assumptions: 


1) om (} 

ii) We neglect the displacement currents in both media, i.e. set 
€, = €, — 0. This means that the transient response is valid for 
times f¢ € <0. 

iii) The dipole is situated at the origin of coordinates and the induced 
field observed in the plane z 0. 


Substituting (19) in (16) and (17), and inverting with respect to the 
Hankel transform, we find 


| MTs) [ o0/72e—@ 
j2 = e } dé . 
H’ se | o 5)0(e) 10, (20) 
: UT(s) F o0(%2—9 
62) JS,(@p) dé, 2 
Nip to | oe i} YP) (21) 


vyhere now 
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Evaluating the above integrals by Gordon’s method and inverting the 
resulting expressions with respect to the Laplace transform, 





; J — 2 —e—p 
H',. 1 G | Ts) l 1 : d (’ é | Ie ds. (22) 
8a* . p b*pdp\ p , 
rf 1x : 
, ] F 7h l > 
H,, soa 3 ("5 an | T'(s)I,(as') K,(as')e* ds. (23) 
f p\p dp . 
where G ‘2 (p : } - = o8, 4a* = op’. 
p dp dp 


4.1. Application A 

A magnetic dipole acting asa step function current source over a homogeneous 
earth 

In this particular case 7'(t) H(t), the Heaviside unit function. Applying 


0 — 





08 


0-6 


0-4 











01 02 03 04 05 O06 


a Laplace transform in ¢ we find 7's) s-!, Substituting for 7'(s) in 
22) and evaluating the complex integral by tables (4) we have 
| (24) 


») >) 1 
Hi, - M HC = 2 d (: oP ‘eal a (“)?. 2/4 
‘ 4a | p op® ap dp p 2 2th 7 | 


where a2 op”. Performing the differentiations we find 
p g 





Hy, = 


A(t), (25) 
4irp* ( 








where 


The fi 
evalué 


we ha 
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= [1—19f an) * [fe 


d 
ert a ert x. 
dx 


where 


[he function A(t) is shown in Fig. 1. Substituting for 7'(s) in (23) and 
evaluating the integral using the Faltung theorem for Laplace transforms 
we have t 


M { « ( * e—x*/87 = 
H M | f- l | | 12 ) dr. (26) 
d tor\dp pdp | T 


. 


4.2. Application B 


A maqnetu dipole acting as an impulsive current source over a homogeneous 





uth 
In this case 7'(t) = S(t). Applying a Laplace transform we have 
T'(s) 1. 
10 —E SS 








The expressions for the induced fields are, after evaluating the complex 


} ri 
n,, — 2! a(}(4 ‘ert( 3.) (27) 


2710 \p dp p 
H', al Sud (dh | POT bt (28) 
trt\dp\p dp} | 8 


integrals 
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Performing the differentiations we have 
, M ; M 
Hi. =-— - Bit), Hi, 


27ap° 270p 


B(t) Hert| ") | 9( lh “( ") + a, ") 
2! oft) Tle of! 


C(t) 16z[(22?+-2z)J,(z) —2(22-+-2+-1)1,(z)Je~, 
- x 
St 


C(t), (29) 


where 





The functions B(t) and C(t) are shown in Figs. 2 and 3. 














4.0 
3.0 
2-0 
A 
c(p| '° 
0 
oF 02 03 04 05 
——> 4. 


5. The induced magnetic field over a two-layer earth when the 
conductivities are nearly equal 
We now make the following assumptions: 
(i) o, = 0. 
(ii) We neglect the displacement currents, i.e. set €, = €. = €, = 0. 
(iii) The dipole is situated at the origin of coordinates and the induced 
field measured in the plane z = 0. 
(iv) The conductivities of the two layers are nearly equal, i.e. 
Ac = o,—03 a, 


») 


Inverting (16) and (17) with respect to the Hankel transform and 
following the method of (3) we see that the resulting expressions can be 
reduced with sufficient accuracy to 
ry* : TI’ ry" 

Ay. = 4,443. 
ie -_ 
Hy, = H 


1p 


+ Hi. (30) 
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The terms H\. and H,, have been evaluated in section 4 and the terms 
arising from the inhomogeneity in the earth’s conductivity are 


» _ MBs) fF ol%—1: 
A". (s) | g2( 72 "3 2472,J,(Ap) dd, (31) 
a No 3 
1} aa 8 re ° : 
Ar MTs) | @°(" m3 2inz J, (Op) dd. (32) 
477 p H27 Ns 


\s in (3) we can write with sufficient accuracy the above in the form 


| AcMT(s)s ~ e720 
i". ee | e-_J,(Op) dé, (33) 
l67 : n* 
0 
» — AoMT(s)s [ype 
H ; = 2 eg" J, (Bp) dd, (34) 
O77 . 
0 
where 7 }(2+ 3) — (6°+-es)', and o denotes the average conductivity. 


Inverting (33) and (34) with respect to the Laplace transform and 


reversing the order of integration we have 


oM 
H’ AoM | 62.J,(Ap) (8, t) dé, (35) 
l6zo0 
0 
oM f¢ 
H; =e . 62), (0p) D(0, t) dd, (36) 
' O70 
0 


where 


(6,1) [ ~@* expl—aio+p)!-+st] ds, (37) 


s+ 
A? tho, B 6? Jo. 
5.1. Application C 
A magnetic dipole acting as a step function source over a two-layer earth 


As before we take 7'(t H(t) and find 7'(s) s-l, Substituting for 


T(s) in (37) we have 





20\1 
(0, t) f H erfo(" 28. (38) 
t / 
Therefore (35) and (36) become 
P al . : 1 r 
H' = rerfe(™ 2) Jy(Op)e® dd, (39) 
O7T0 
0 
ol . [h?o\} r 
a. on *erfe(“ ~e | 0J,(Op)e-P* dé, (40) 
f l670 dp. 
0 
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where 7 . (Pa) By (5, p. 394) 
pdp\' dp} * 
\1 fp 2 1 2) 
Hi}. - —= *erfo( “2 “DI mie ’ (41) 
32 \zot | t St] 











i A a 


www iu 
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and by (5, p. 393) 


m h . fh2o\1 : 
H! _ Ao L fel” o\hd | x/a0 (49) 
' 327 dp 
Performing the differentiations we have 
” M ” J y 9 
H;. SS re __ Ae ! FW), (43) 
. 4rrap* ? l6zp 


where 


The fu 


5.2. A 
A ma 


In t] 


36) gi 


H;. 


H;, 
Eva 
H;. 
H;, 
where 
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ynere 


D(t) erfe( 7°) * (2z—1)J,(z)— 221, (z)](272*)-*e~*, 


re) = er) (Ge 


op” e > ° 
Kh? p. 


bal / : 
The functions D(t) and F(t) are shown in Figs. 4 and 5. 


5.2. Application D 


{ magnetic dipole acting as an impulsive current source over a two-layer earth 





In this case 7'(s) |. Evaluating the integral (37) we find that (35) and 
16) give ( 
\oM h?o\4_ h{(o\3 ., [h?a\h P 
H; lence / erfe| }" J i;) | ; "| J, (Ope Bid@, (44) 
0 
\oM d[_, [h?a\}.. hfo\3__.,[h2o\4] 
; rfe <¢ “erf’ . AJ,(Op)e-Pt dé. (45 
Hy | 6770* ¢ at ss | t : (;) - ( t | ol py — 
6 
Evaluating the integrals as before we have 
ol = h2o\} h(o\2 ,,fh?o\3 ; 2 
H; : ‘( . }* 7) erfe| “)*9 (°\? ert (" “}* (2 , (46) 
32 \no%t) ~ |. t 2\t z Sf 
oM d . [h?o\4 h{o\2 ..,{h?o\} ” 
H, wines erfe| i (7) er! (“ > at re, (47) 
j 3270 dp \ t ? t t 


where iy ap". 
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THE DIFFRACTION OF A DIPOLE FIELD BY A 
HALF-PLANE 


By BETTY D. WOODS 
(Department of Mathematics, Manchester University) 
[Received 7 March 1956] 
SUMMARY 

The problem of the diffraction of an electromagnetic dipole field by a perfectl 
conducting half-plane is discussed. Solutions for arbitrary orientations of the dipok 
are determined by extending a method due to Bromwich, which was based on thy 
solution of a scalar problem and valid only when the axis of the dipole is paralk 
to the edge of the screen. Results are given for an electric dipole with its axis 
normal to the screen and for an electric dipole with its axis normal to the edge of 
the screen and lying in a plane parallel to the screen. 

1. Introduction 

THE problem discussed in this paper is that of the diffraction of the electro- 
magnetic field of an oscillating dipole by an infinitely thin, perfectly 
conducting half-plane. This problem has been discussed by Senior (1), 
who obtains a solution by expanding the dipole field in plane waves. The 
method can be used for any position and orientation of the dipole but the 
calculation is laborious. Senior gives an explicit solution only for an 
electric dipole with its axis perpendicular to the plane. Another method 
has recently been given in a paper by Heins (2).+ 

Bromwich (3) has given a simpler solution for a dipole parallel to the 
edge of the conducting plane based on the solution of the corresponding 
scalar problem, that of the diffraction of the sound waves of a point source 
by arigid screen. The purpose of the present paper is to extend Bromwich’s 
method to arbitrary orientations of the dipole. 

The solution of the scalar problem is the determination of a function ¢ 
which satisfies the wave equation everywhere except at points on the screen 
and at the singularity which represents the source. Also it contains only 
diverging waves at infinity and its normal derivative is zero on the screen. 
There is a similar problem in which 4, instead of its normal derivative, is 
to be zero on the screen. Solutions to both problems, which are closely 
related, have been giver: by Macdonald (4), inter alia. 

Bromwich’s method is to take this function ¢ multiplied by a constant 
unit vector in the direction of the axis of the dipole as the Hertz potential 
for the electromagnetic field problem. The resultant field satisfies Maxwell's 
equations, has the correct singularity at the position of the dipole and 

+ Note added in proof. Investigation of this problem has also been carried out by Yu. 
V. Vandakurov, Zh. Eksp. Teor. Fiz. 26 (1954), 3-18. 
(Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 1 (1957)] 
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contains only diverging waves. By using the value of ¢ which satisfies the 
sppropriate boundary condition for each of three mutually perpendicular 
orientations of the dipole the tangential component of the electric field on 
the sereen can be made to vanish. Thus Bromwich’s method leads to a 
field, for arbitrary orientation of the dipole, which satisfies Maxwell’s 
equations and the appropriate boundary conditions. 

[It is known, however, that for diffracting bodies with sharp edges these 
boundary conditions are not sufficient to determine a unique solution of 
the field equations. Recent investigations by several authors (Copson (5), 
Meixner (6), Jones (7), Heins and Silver (8)) agree that the correct solution 
to the physical problem is defined by the order of the singularities of the 
various field components near the edge. It is only for an electric dipole 
with its axis parallel to the edge of the diffracting half-plane that Brom- 
wich’s solution satisfies the edge conditions. 

The difference between any two solutions of the field equations which 
satisfy the boundary conditions on the screen and at the source is a field 
which satisfies Maxwell’s equations everywhere except on the screen and 
satisfies the boundary conditions there. The method adopted in this paper 
is to determine the solution of this “~homogeneous’ problem which, when 
added to Bromwich’s solution of the diffraction problem gives a field 
which satisfies the edge conditions and which is therefore the correct 
solution of the diffraction problem. In this way the problem involving 
either an electric, magnetic, or sound dipole can be solved. The details 
are given here for two cases, one being that of an electric dipole with its 
axis normal to the screen, for which the solution is shown to be the same 
as that given by Senior, and the other that of an electric dipole with its 
axis normal to the edge of the screen and lying in a plane parallel to the 
screen. The solutions for these two positions of the dipole, together with 
that given by Bromwich for a dipole parallel to the edge of the screen, 
can then be used to determine the field of a dipole orientated in any way. 

I wish to express my thanks to Mr. D.S. Jones for suggesting this problem 
and to him and Mr. E. Wild for their valuable advice and criticism. 


2, Bromwich’s method 

Suppose that an infinitely thin, perfectly conducting half-plane occupies 
the region y 0,2 > 0 of a cartesian coordinate system O(x, y,z) and that 
either an electric or magnetic dipole is situated at some point in space, its 
direction being represented by the unit vector u. 

The electromagnetic field must be such that: 

(i) It satisfies Maxwell’s equations everywhere except on the screen, and 


it the dipole where there must be a dipole singularity. 
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(ii) It represents diverging waves at large distances from the origin. 

(iii) The tangential electric field component vanishes on the screen. 

(iv) In terms of cylindrical polar coordinates r, 0, z, the x and y com. 
ponents of the electric and magnetic fields are of order r- and the 2. 


1 


components are of order r} near the edge of the screen. 
For the case of an electric dipole Bromwich’s solution 
EB) (x, Y, ze ickt HP®) (x, y, z)etckt 
is given in m.k.s. units by 
&4 = (VV.+k?)du, A) — ikY(V A bu), 
where Y is the intrinsic admittance of free space. 

The function ¢ is the solution of the wave equation (V?2+ k?)¢ 0, which 
either vanishes or has vanishing normal derivatives on the screen and js 
such that near the position of the dipole 6 ~ e-*®/R, where R is the 
distance from the dipole. 

€ and A”) satisfy (i) and (ii) and can be made to satisfy (iii) by using 
the ¢ which satisfies the appropriate boundary conditions for each position 
of the dipole. Thus for the dipole parallel to the x-axis 

a9 ay J 
: 7d = é C2 
G(B) : ‘: kd, GP) $ : 
Ox" ‘ CXLOZ 
6%) — €6®%) — 0 ony = 0,x > 0, when d = 0 there. 
For the dipole parallel to the y-axis 
- i 
aiB) — ¢ ap) — “ ? 
; Cxoy “ eyez 
(B) (B) _ . > io ¢ ‘ ——-, 
GC é\ 0 ony = 0, x > 0, when éd/éy = 0 there, 
and for the dipole parallel to the z-axis 
a9 a9 
od ;, OS. x. 
GB) — —? E(B) P41 kd, 
OXCZ oz* 
é%) — 6% —) ony= 0,2 > 0, when d = 0 there. 
For the case of the magnetic dipole Bromwich’s solution is given by 
ik 


Ee) p(VAgu), A — (VV. +h*)pu, 


where again €” and A”) satisfy (i) and (ii) and satisfy (iii) when ¢ is } 


suitably chosen for each position of the dipole. 
For the dipole parallel to the x-axis, 
EB) — 0) E(B) th Op 
7 ° Y ey 


é) = 0 ony =0,x > 0, when éd/éy = 0 there. 
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n For the dipole parallel to the y-axis 
Or 6B " ad 6) a ¢ 
; Y oz 7 Y ox 
EU E(B 0 ony = 0,x2 > 0 when ¢ = 0 there, 


nd for the dipole parallel to the z-axis 


— 
y ] vk . GAB) 0 
Y oY 7 
GU 0 ony 0.x 0 when 6¢/éy 0 there. 


3, Behaviour of Bromwich’s solution near the edge of the screen 
We consider the case of an electric dipole with axis parallel to the y-axis 


und situated at the point (ay, ¥»,2)), for which Bromwich’s solution is 


sl 84 td a2 
; am — (2? CP, pag © 4) (3.1) 
\¢ voy cy~ Oyoz 
KB) — iky =, 0,2), (3.2) 
Oz Cx 


here d = ¢,+¢, and in the form given by Senior (1), 
Ms 


: | H'?(kS cosh p) du, 


f] 


HY (kR cosh p) dp, db 


where R p27? 2rr, cos(8@—O,)+(z—zZ,)?|? is the distance from the 


) 


lipole, and S r?+-r2— 2rr, cos(O+-6,) + (z—Z»)?|! is the distance from the 








image of the dipole in the plane y = 0; ry, 6,2 are the cylindrical polar 
coordinates which specify the position of the dipole. Also 


(# H ) 
0 > 
cos : Ks sinh 


1 2, (797) (6 T 69) 
R 2 : S ) 


COs 


The solution of the homogeneous problem is most easily determined by 
ising Fourier transforms, and therefore to determine the difference between 


Bromwich’s solution and the correct solution we introduce the Fourier 
transforms E™), H) with respect to z of &, A defined by 


~ F ; l i ; 
Ee , EWe-is= dz, A H Boise dz, 






nd we consider the behaviour of E“” and H™ near the edge of the screen. 
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It is seen from (3.1) and (3.2) that 


Ew — (29 &® | p2@ 5.2% (3.3 
Gxoy’? ey? by | . 
H® — inY -is®, 0, ae) (3.4 
ox 
l f ia 
where do — | de-'** dz. 
y(27) , 


Now © is easily determined, for as ¢ satisfies 
V7*b-+ k?h+-8(a—ay)5(y — Yy)8(z—2) = 9, 


® must satisfy 


: ; ¢ iszo " fe a 
V3O+ 2+ 8(e—ar9)5(y—Yo) = 0. (3.5) 
\ (27) 
ae @ 
where V3 —.+—) and 6 is the Dirac delta-function; also «? = k?—s' 
Ox* oy" 
and we define « as the branch of (k2—s?)! for which « > 0 for |s| < k and 
ix > 0 for |s| > k. Equation (3.5) is the same as the equation satisfied 


by the magnetic force, e'**,/(27)®, of a line source on the line ry, 4) in the 
. N 0 0 
presence of the screen and from the solution of this problem given by 
MacDonald (4) we find that 


© = 0,+9,, 
where Er 
ep —t8zo a es 
ri) e—ikR, cosh p du, (3.6) 
1 } (27) 
vA ms e 
Hs, 
e820 » es - 
o, p—ikS, cosh p dp, (3.7) 
- I (27) 
vy aml e 
R 72+ r2— 2rr, cos(@—8@4,) |}, S, = [r?+r2—2rr, cos(6+ 4), 
1 0 0 0 1 ™ . 
1 a f2v(ror) (0-8 mht 2V(%0") cos (4-50) 
Si sinh- V\'0 cos : o) ; Re. sinh-! v("o cos 5 Ory. 
1 = ae ‘ 


Equations (3.6) and (3.7) can be verified from the inverse Fourier transform 
using the substitution 

t = —tx(r+19)+28(z—Zp) 
given at the end of section 4. 


_— e ° — ° . 1(B) 
lo derive the edge singularities in the transformed components EY”, 
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B®), EB, and H®), H\»), HY) it is sufficient to determine the leading terms 


nthe expansion of ® in powers of r! near r = 0. Thus 
0 
pe —t8Z9 wy e—'tKro—iszo Pr 1 fal -@ 
0, = | é KT COSD ft du — a ( J eos: To) O(r), 
ky (<7) , k an) = 
0 
? 8 ° ikro—iscto Ir 1 0 l g 
®, e—ixrocoshw dy 4 . — ( JFeos' . 0) O(r) 
k,/(27) , k TI 2 
is -0. Hence 
0 
, 8 2\1 e 2 ikrg—iscze Py | 1 fa fa 
ty) }* | aa Daas du be ( ; }* cos ®cos—+ O(r) 
k T k \71'o) 2 2 


Using this result and carrying out the differentiation in (3.3) and (3.4) we 


find as r > 0, 
0 
. CU 36 | , 2\3 ¢ 
k DB = sin ~-+ O| } H\? s) oe —'820 J: | ¢ ixrgcosh pu du O(r+), 
1] C 36 | B 
Ri cos { O| ) AY” 0). 
2) 2 ? 

- isC ikY C 0 
he sin 6+ O(1), H®) cos = + O(1), 

Yr 2 


4. The homogeneous problem 
To eliminate the incorrect edge singularities in Bromwich’s solution we 
ave to find an electromagnetic field €@e%e*, AOeick which satisfies 
Maxwell’s equations everywhere except on the screen, represents only 
liverging waves at infinity and has a vanishing tangential electric field 
omponent on the screen. We consider only solutions & and A© which 


have Fourier transforms E® and H with respect to z. 


’ 


The equations to be satisfied by the Fourier transforms ZO, BX), B&, 


nd HW), HH’, H© of the field components are 

rk ‘ C Hy) \¢ 7 ° ry rmwe (C) ° ’ 
a Ee” iskt KY RO — fe” gq, 
} oy oy 

iL ( a Rt of (C) 

: H is E\ at ’ tkY EW) isH© H, , 
} Ox p Cx 

Lt 3B aR i (C) aH ©) 

kH K\ Eg py po — Cf on; 
} Ox cy ” Ox oy 
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From these it can be seen that all the transformed components can be 
expressed in terms of E&, H% and their derivatives. Thus 


RO isOE™ wk eh HO iYkoE® is eH©® 

ps = s > : 3 r = - t > ; 

Pi x2 Ox Yu? oy ; K2 Gy K2 Ox 

RO) isOE® § ik dH® He - Yk oR is dH) 
yu Ke ey Yu? Ox * Kx Ke ey 


(4.] 


It can be seen from (3.3), (3.4), (3.5) that BY, BY), HY, and H\) can bk 
expressed in the same way in terms of EY and HY). From (4.1) and the 


boundary conditions ZO) = ES’ = 0 on the screen there follow th 
boundary conditions 
“cn _. OH? 
i = —* 0 on the screen. 
oy 


As & must satisfy (V2+ k2)66 = 0then EO must satisfy (V+ x?) BO = 0) 
and the only solutions of this equation which satisfy the boundary con 
dition and the condition that ES ~ Qe-‘«"/r! for large 7, Q a constant 
(which is necessary in order that the solution of the original problem should 
contain only outward travelling waves at infinity) are 


ners “ as . no 
BO) = > A,,, H®(«r)sin = (4.2 
n 


where v has integral values and the coefficients A,, are arbitrary functions 
of s. We define H?)(«r) for « real or imaginary as a Bessel function of order: 
which has the asymptotic form (2/z«r)e-*"-!"7—-!™ for large r. 
In the same way for H“” we obtain 
H® = Y B,, H®(kr)cos . (4.3 
— ; 2 
where the coefficients B,, are arbitrary functions of s. 

We consider now the electromagnetic field whose Fourier transform is 
given by E = E”+E©, H = H®-+-H, and we determine A,,, and B,, 
by equating to zero the terms of order r~! in the expansion near r = 
of E, and H,, which are thus made to satisfy the edge condition. It follows 
from the relations (4.1) for the B and C fields that the other components 
of E and H will also satisfy the edge condition. 

For the values of the coefficients we find 


me | = 
A, (SY Se ata eve)’. 
! 5 ! > 
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Hence | 
k 17? 
Hi! Ye-ixtr-+re ae 
7 
and therefore 
/ CD ise tk +r —( 
y is 
cy k oi 
sigs 
H ikY < 1Ye—tKer 
Cx 


From the definitions (3.6) and (3.7) and t 


a) and (y—y,) and S, a function of (x 
f sily that 
c® C (®, QD, ) gtk +1o)—t829 
cy CY ko 
Pali) c® e—tKkr+rq)—t829 [ 9 
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) 1 

cos ban? 
Pu 2 2 


2 G : 
cos sin —, 
5) 9 


» ] 
e)—iae - by) Gy _o 
oe tae cos — COs —. 
TY 1 2 Z 


he fact that R, is a function of 
%») and (y+y,) it follows 


» 1 
( 2\; 6. 6 
- - cos sin-—, 
ST >) 9 
TH 7 zZ ya 


2 
“COs cos —. 
5) 9 


Using these relations we can express EH, and H, in the form 


1s (, 
c Yo 


®,), 


For the other components of E and H 


H, 





,- od 


OX, 


iky 


we have, using (4.1), 


; ike ix(r+ro)—iszo 2 1 fa . 
EB (®,—,) mee | “ cos — sin 
OXLCY, od K 719 r 9 2 
’ 0? ‘ ike—tur+re)—tsz0/ 2 \} 6 ] 
c= - (®,—®,)+ k?@ — —— . “cos—cos-, 
CYCYy K Ty? 2 2 
’ af F K(?'+19)—i8Z9 9 1 6 ) 
H skY@O-+ J = ( “ cos—2cos—, 
COX K TY 1 2 2 
gic 2 \4 G « 6 
i =] | “cos sin -—. 
C2 K | \7r7' 91 2 2 


Finally, by taking the inverse Fourier transforms of E and H we obtain 


the correct solution €, A of the problem. Thus for the x-component of 


the electric force &é. we have 





ds 


K 


e—ix(r+ro)+ is(z—Z9) 


(4.4) 
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To evaluate the integral in (4.4) we make the substitution 


t = —ix(r+75)+is(z—z,). 
Then 8 = —it(z—2)+(r+1r9)(k?U?2+#?)!, 
where U = [(r+19)?+ (z—29)?}}. 
e! dt 


rhe integral becomes UL 
a 
where C is a contour which can be deformed into the imaginary axis taken 


from —ioo to ico with (k?U?+#)! > 0 for |t) << kU and i(k2U2+#)! > 6 


for |t} > kU. This contour can be further deformed into a Hankel loop 
round the branch point ¢ = —ikU, which gives for the value of the integral 
7H (kU). Therefore 
oc ik - . ee 
é, = —— (¢.—¢,)+- H® (kU )cos sin —. 
Ox0y,  ~ (797 2 2 


The other field components can be determined in a similar way. Hence 
for the complete solution we have 





= “7, } 
6, = (bby) +, HKU eos sin, 
‘ ‘ 2 1 ee 0 5) 9 
CXCY,) (77) = a 
- “I. J j 
6, = (b—$)) +h ——"* HEU) c08% eos, 
CYCY, V (797) “ . 
6, - it ( 2 —d¢y), 
CZ0Yo 
#, = ay ®t et kY (z—29) H®(kU)cos 4 wien, 
é Cc Zo i \ (775 r) 2 2 
.) — ) 
Free kY (z 2) H® (kl T)oos “sin : 
y l (fof) 2 2 
\ (79 
#, = —ikY Op 
Xo 


It is seen that €, # in the above form is the solution given by Senior for 
this problem. 


5. The dipole with its axis parallel to the x-axis 

For the case of an electric dipole with axis parallel to the a-axis and 
situated at the point (2, ¥,2)) the procedure is the same. Bromwich’s 
solution is 


gm — (2%, p24, 28 =r) HED) ity (0, &@, |, 


C2 CxCYy CXOCZ oz oy 


where now ¢ = ¢,—4¢y, 


? 
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The Fourier transforms of € and A“) are 


oD . eo . e® —_ , ed 
Eu ae ek (0. we, —%), 
Ox* CxXCY OX | oy 
vhere @ ®, ®,. 
Near the edge of the sheet 
' isD . 6 me ikYD 6 . 
Ss sin + O(1), H‘®) —eos=+O(1), (5.1) 
r 2 : r? 2 
ker e-txre isco 2\h . 8 
where D | )sin a 
. ana ») 
k 119 2 
The coefficients A,,, and B,,, involved in the homogeneous problem are 


letermined by 
l€ scribed. 


Thus we find 





Hence 
isSé 
BR 
H Ye 
en for E. and H. we obtain 
: . e® LS¢ 
k 1s 
Cd 
malt) 
Cu 
} which can be written as 
ID) is (D, 
OXy - 
By using the relations (4.1) 
found. Thus 
E —— 5 
ue Lo a , 
> ID (®,—9,) 
| {CXo ™ 
| r ¢ m 
H } 
} Crg K 
| 





eliminating the terms of order r-! 


in (5.1) in the way 


a 
an. 
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‘ 1 
2\i. 8, 
sin 
TV)? 2 
2% 
To ~~ a 
TY )1 


8 0 ’ 
sin cos —. 
» » 


mrivd—te/ 2 \4 . O, . GO 
; “sin sin -—, 
k 7 1 2 2 


0 
COS —, 
9 9)’ 


; 2\34. 6 
iY e—tx(r+ro) fon J*sin 
79 r 


®,), H,=iky(@, 


CYo 


{ ®,). 


the other components of E and H can be 


k?@ — : sin — sin -, 
K 7119 2 2 
ike-ixtr+re)—iseo/ 2 \2 . Q 
“sin COs -—, 
K 71 91 2 2 
2\4. @ 0 
| }*sin COS —, 
71 2 2 
ix(r+1ro)—isz 


a 
si sin —. 
9) 


1]. 86 
fxn : 
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Hence by taking the inverse Fourier transforms of the components of f 


and H we have for the solution 


é ba. + k*d- os H®)(kU)sin ae . 
. Oxda, | V(rer) ° ; 2 

é, te 4 H®(kU’)sin _— 4 
CYyOX, (77) 2 2 

bn ee 
. 0z0Xy 

FF . kY (2 ~Zo) H?(kU)sin econ 4 
: Uy (7) 2 2 

H, = —izy &. aes =a) HEU sin sin ®, 

CZ J (To? 2 2 


#,- ar ar. 


CYy CYo 


The field due to a dipole orientated in any way can now easily be deter. 
mined from the solutions for the above two cases together with Bromwich’s 


solution for an electric dipole parallel to the z-axis. 
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THE PERIODIC SOLUTIONS OF THE DIFFERENTIAL 
KQUATION OF A RESISTANCE-CAPACITANCE 
OSCILLATOR 


By A. W. GILLIES (Northampton Polytechnic, London) 
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SUMMARY 
\ third-order differential equation is considered of a form which arises in con- 
exion with a resistance-capacity oscillator, the equation being normalized to the 
V6 | “ 
D 1)(D - } eD ja g(D) > Cap” Lyn g(D)2Bcoswt, 


vo n 


where € and yz are small parameters which in the main part of the discussion are re- 


ited by « u?, and g(D) is a particular polynomial operator of the third degree. 
[he procedure previously applied to a second-order equation with unsymmetrical 
n-linear damping, is used to obtain the periodic solutions having the period of the 

forcing term when w is near to 1, i.e. for the case of fundamental resonance. It is 

hown that the resonance curves are of the same form as those obtained for the 
md-order equation and that the stability of the periodic solutions is determined 
variational equations which are likewise identical in form to those previously 
tained for tl second-order equation. 


1. Derivation of the differential equation 
THE equations of a triode oscillator may be written in the form (cf. (1)) 


v e,+aZt,,, 


] 7 


Zi, é v 


a—Va 
where v,, v, are the anode and grid voltages, Z is the impedance of the 
inear circuit at the anode terminals when the grid terminals are open, « is 
the voltage transfer ratio from anode to grid through the linear circuit, 
ind e,, ¢, are the effective e.m.f.s in the anode and grid meshes of the 
inear circuit. 

The currents and voltages are increments from the quiescent point, and 
t is assumed that the circuit is biased so that no grid current flows. 

If it is assumed further that 7, is a function of mv,+v,, with m the 


( a’ 


implification factor of the valve, and if we set 


? mv, + Var 
f me, + Ca 
then we obtain (l—ma)Zig e—v. 


[Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 1 (1957)] 
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If we refer to the series C, shunt R circuit, Fig. 1, we obtain 


a Ro{(m+1)R°46R?°X+5RX24 X93) 
RGR? +4 RX +X} +1 R84 6 RX +5 RX? NS) 








with 1/X = jwC, from which we get 
Ro{(m+ 1) R3+6R2X +5RX24 X3}7, 
= {R,(3.R?+4RX + X?)+ (R3+6R2X +5RX2+ X3)}(e—y), 
Since we have been concerned only with the linear circuit we may obtain 
the differential equation by replacing R/X = jwRC by CRd/dt, and by 

















Fic. 1 


a change of time scale this may be replaced by d/dt = D. If further R, 2 
is assumed negligible in comparison with unity we obtain the differential 
equation 


{(m+ 1)D®+6D?+-5D-+-1}R,i, = {D®+6D?+5D-+ 1l}(e—»). 


j 
in which e will be a known function of the time (in particular a simple 
periodic function) and the valve characteristic determines a non-linear 
functional relation between i, and v. It is usual to express 7, as a power 
series in v, but in the present case the number of parameters in the final 
equation is smaller if v is expressed in powers of i,, and this will be done. 
although the difference is not important. It is assumed that such a repre- 
sentation is valid over the operating range of the characteristic. 
Write x, = R,i, and y instead of v, so that we have 


{(m+-1)D®+6D?+5D-+ l}a,+{D®+6D?+5D-+ ly 
= {D®+6D*?+5D-+ le. 


Suppose that in the operating range y can be expressed as a series in 
powers of x,, of the form 


_ - to’ 22 ’ 273] 
Y = C1 XCo PU Cg bX + 


in which ¢, is positive, and in a normal valve characteristic c will also be 
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positive, » is a small parameter characterizing the magnitude of the non- 
linear terms. The equation then becomes 
m ; : , y ; eo " c n— lyn 
(5 1| D? +6 5D+ 1}, -{D?+4+6D?+5D+1} > <n 
Cc | 
tl 1 


1+-c, 


(D3 +6D2+5D+ 1}. 
Ite, 


If the non-linear terms are omitted (u = 0), and we set e = 0, we are 


eft with the linear equation 


| Vt mee 3 , 
(; i 1|D*+6D® 15D) at 0. 


l 
This equation has a simple periodic normal mode if 


m 29 
l+c, 


when it reduces to (6D?+-1)(5D+1)a, = 0. 


This represents the condition of critical regeneration. The linear equation 
has an exponentially increasing normal mode ifm > 29(1+c¢,), a decreasing 


mode if the inequality is reversed. Set 


Wi 


und 302, x: 


the full equation then becomes 


; rn bs ~ c n-1 xa \" 
((1-+e)D3+1D?+4D+d}e+{D3+6D?+5D4+1 5 2 ; 
; — 1+¢, 30 
n=2 


(D3+6D?+5D41 
+e, 


\ further change of time scale, replacing dt by v6dt, brings the equation 
to the form 


| 6 boon _ ¢ n-1 a\” 
|(D2-+-1){ D+ ~ <D\x-+{D?+6v6D2+ 30D +66} > — A fo 
| | 5 1+c, \30 


n=2 


{D8 +-6,y6.D?+ 30D+6v6} 


t 
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It is convenient to associate a numerical factor with the operator 
{D*+6V6D?+ 30D-+ 6v6} and set 


299(D) = D®+-6v6D?+ 30D+6Vv6, (1) 
20°" 
: 29) ™ t; == 2a 

30"(1-+c,) 
, : 29e 

If we further assume iw = 2Beos wt, 
1¢ 
1 


the equation becomes 
6 
(D2 (D+ =) DY +9(D) ¥c, pa" = g(D)2Beosut. (2) 
0 n=2 


The operator g(D) on the right-hand side of this equation produces a 
variation in amplitude and phase which will be small if only a small fre- 
quency interval is considered. This factor may therefore be absorbed by 
replacing g(D)2Bcoswt by 2B, cos wt. 

The equation is then of the general form 


{(D?+1)(D+a)—ek(D)}a+-9(D) ¥ c, w"—1a" = 2B, cos wt, (3) 
n=2 


in which a is a positive constant, « and jy are small parameters, and k(D), 
g(D) are polynomials in D, of degree not exceeding 3. 

The equation is thus of the third order, containing small non-linear terms 
dependent on the parameter j., and such that if these, together with small 
linear terms dependent on «, are omitted we are left with a linear equation 
with constant coefficients having a simple periodic normal mode of angular 
frequency normalized to unity together with a decreasing exponential mode 
of time constant 1/«. 

The present paper will treat the special equation (2), though it will be 
clear that the method is applicable to equation (3). The generalization to 
a wider class of equations which will include (3) is however deferred to 
another paper. 


2. Outline of the method 

Weare concerned with (2) where g(D) is defined by (1). Thec,, are assumed 
to be O(1); « and yp are small parameters which will at first be treated as 
independent, but later they will be related by « = p?. In the case worked 
out in detail B will be taken as O(e) and (w— 1) likewise, i.e. the amplitude 
of the forcing term is small, and its frequency near to resonance with the 
simple periodic normal mode of the approximate linear equation. 

In the equations of the first approximation only c, and c, occur, so that 
effectively the valve characteristic is approximated by a third-degree poly- 
nomial. 
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The method of solution will follow closely that used in connexion with 
the equation v—(e- pkv p*v?)b +o = 2Bw, sin w, ft, 
which was solved in a previous paper (2). The procedure for determining 
periodic solutions having the frequency of the forcing term is exactly the 
same; but in considering stability account has to be taken of the fact that 
the present equation is of the third order. 

In a normal characteristic c, is positive, and if € is positive, regeneration 
s beyond critical. It will be shown in this case that, if the ratio c3/c, is less 
than a certain value, the resonance curves are of the same general form as 
for the second-order equation, and that the stability of the periodic solu- 
tions is likewise determined by variational equations of the same form as 
those previously obtained. If, however, c3/c, exceeds this value, the ampli- 
tude of oscillation increases and the final state cannot be determined from 
the equations of the first approximation. 

If ¢ is negative the regeneration is less than critical. In this case, in the 
absence of the forcing term a small oscillation will not grow from zero. 
With the forcing term, when c3/c, is less than a certain value, the resonance 
urves resemble those of a system with a non-linear restoring force, such as 
acompound pendulum. If, on the other hand, c3/c, exceeds a certain value 
tappears that if an oscillation of sufficient amplitude can be started it will 
grow, and again the final state cannot be determined from the equations 
of the first approximation. 


3. Periodic solution developed in powers of B 


Equation (2) may be rearranged in the form 


2Bcoswt = ((D)x+ Fc, pw" 12", (4) 


| -| cD, 


with C(D) fli/g D)}\(D* (D4 : 
which is of exactly the same form as (7) of (2). We may therefore proceed 
sin (2), writing 
x CY B+ 7 B2+ tam Br... (5) 

Substitute into (4) and equate coefficients of like powers of B; we thereby 
btain the 2” in succession as periodic functions of t of period 27/w. 

Using exponential notation, and reverting to 7 instead of j for ,/(—1) we 
thus obtain 

exp iwt + exp(—iwt) C(D)x® 


4 (2) . «1)* 
O = C(D)x-+-cy pa (6) 
0 C(D)ax®) | 2c, pxa'?) LC, pall? 
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from which aD — gf) 4 ofl) | 
] : l ‘ 
where 2) — —expiut, 22) — ——exp(—twt). 
G1 G4 
Similarly 2) — 2) 4 72) 1 9(2), 
where = - 
22) — - f 2et1* — — “ 2 exp 2iwt 
le Cite 
>) ») 
y(2) ON? (1) »(1) _ _ She, (7 
“— "ies iil, _ y 
So £1014 
Le. uc : 
z), - — Byfty - = 2 exp(—2iwt) 
6-2 6—1$-2 
Again 2) = gS), + o((8) 4 gS), + o((8), 
where 1/202 
— ial 2,.(1)3 
~~ 2 — Cg poe, 
Cs\ Co 
8) 
Ll (2c, . 4c 
.(3) at 2 § 2,(1)? »(1) 
uy eo - | + s 3c} vy vn 1 
>1 So So 


and x), and x) the corresponding conjugates. In these equations C, 
denotes f(-+-riw). 

Proceeding in this way the 2” are determined in succession in the general 
form am) = > (9) 
with 7 running through alternate integers from —n to +n. 2”, is the 
complex conjugate of x!” and the two together give a real term of rth 
harmonic frequency. 

The series (5) then gives a formal solution of (2) which is periodic with 
period 27/w. The series is in powers of B, but, since x contains p”~" as 
a factor, it may equally be regarded as a power series in p. 

An interval of assured convergence is obtained by considering the 


equation ¢ Ny—cp, 7? —cp? 7 —..., (10 


where ¢,| <¢ (es = 8, 3...) 
and N, , are positive constants to be chosen later. Equation (10) with 
N + 0 defines » as a function of €, vanishing at € = 0 and analytic ina 
circle extending to the branch point nearest to the origin given by 

dé 


0, 
dy 


This branch point is at € = R, » = S where 











» ma 
being 


subst 
the 7 


and 


the s 


for a 

Tl 
thus 
and i 
sion 


Tl 


If .V 
hanc 
genc 

If 
rw is 
ditic 
i.e. 1 
cony 


W 


Whe 
dent 
sam 
the 

solu 
deal 
case 


4. I 
\ 


wit! 














ON A RESISTANCE-CAPACITANCE OSCILLATOR 107 


7 may therefore be developed in powers of € provided € < R, € and 7 now 
being assumed real and positive. If this series is obtained by setting 
n= 1ME+ of MEF +-... + fE"+..., (13) 
substituting in (10) and equating coefficients of like powers of € to determine 
the 7”) in succession, it is clear, on comparing with (4) et seq., that, provided 
2B<E&<R (14) 
and NV is chosen so that 
N C(+riw)| (r 0; 4, 2Z....3, (15) 
the series (13) dominates (5) in the sense that 
2m) Bm s b 3 a”) | Be < yfmen (16) 
lor any p f4y4- 

The absolute and uniform convergence with respect to ¢ of the series (5) is 
thus established in a definite interval of sufficiently small values of B and p, 
and it may be rearranged according to frequency to give the Fourier expan- 
sion of the solution. 

The interval of assured convergence is 

i in N\} . 
B IV 4 2e—2e(14 }*). (17) 
244\ ' c | 
If V can be chosen independently of u, and pu, is chosen to make the right- 
hand member of (17) greater than B,, then the interval of assured conver- 
gence is B < B,, for any » < p,; here B, is any convenient constant. 

If, however, w 1. we have |C(iw) e, so that to satisfy (15) when any 

w is near to 1, NV must have a value less than e which is small in the con- 
ditions contemplated. If « is independent of p it is still possible to choose p.,, 
i.e. if the non-linearity is sufficiently small the series in powers of B will still 
converge. 

When JN is small (17) becomes approximately 

B < N?/(4y,¢). 
When NV is O(c) the interval of assured convergence may be fixed indepen- 
dently of the magnitude of the non-linear terms if ¢ is larger than or of the 
same order as y!; but if e is of a higher order of smallness in relation to p, 
the interval of convergence will tend to zero with ». In this last case the 
solution in powers of B fails if any rw is near to 1. This situation may be 
dealt with by the methods of the following sections in which the special 
Case « u2 and w near to | (i.e. fundamental resonance) is treated in detail. 


4. Periodic solution in powers of b 
Whatever the values of ¢ and p, (17) determines a definite interval of 
convergence, provided e is not actually zero, so that for a sufficiently small 
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value of Ba solution of period 27/w certainly exists. (The stability of such 
a solution remains to be discussed.) Let the fundamental mode of such a 


solution be 2b cos(wt +d) 


and replace (5) by 


with a) = aY+a0 = expi(wt+d)+expi(—wt—¢). 


x= ONH+ a@H2+... + alMHr+... (18) 


Following the procedure developed in (2), substitute (18) into (4) and 
choose the x, 2®,... to annul all terms in the successive powers of b of 
other than fundamental frequency. When expressed in terms of 2 and 
2), we obtain the same expressions for the 2” as before, with the omission 
of x” and x", when v is odd, and the omission of terms derived from these 
in the higher orders. Thus the terms of 2 are still given by the first 
expressions in (7) and 2 = 2)+<2®, with 2% given by the first of the 
equations (8). 

The terms which were previously cancelled by the x{") and 2™,, now 


remain to form an equation containing only terms of fundamental fre- 


quency, i.e. containing either expiwt or exp(—iwt). These two sets of 


terms are conjugate complex numbers, and the full equation will be satisfied 
if we satisfy the equation formed with one set of terms only. Those in 
expiwt give the equation 


. sa 2c2 = 42 rakes 
Bexpiwt = C(iw)a\! b- ( F fe ae Beg) wa” be +... (19) 
Ce So 
Since a) — expi(wt+d) 
this may be rewritten 
om acs te P 
Bexp(—id) = C(iw)b- ( wt Sens Beg) 2b 
Go So 
the right side of which is an infinite series of the form 
Zz Az, - perpen +1 
n=0 
where the A,,,,, are functions of the c, (r < nm) and of the 
C(triw) (r | ae ae 3) 


when n > 0, while A, C(tw). 
The relationship between the series for x and 7 remains valid provided 
2b < PE < R/N 
and N satisfies (15) for r = 0, 2, 3,... but no longer for r 1. NN can there- 
fore now be chosen independently of «; and so, given b,, we can always 
ind py, so the > series (18) converges for b < b, provided p < yp. 
find yw, so that the series (18) converges for b < b, provided ; by 
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The series on the right of (19) is likewise dominated by the 7» series 
multiplied by N and will therefore converge under the same conditions. 

The possibility of finding a periodic solution of (2) now depends on 
obtaining real solutions of (19) for 6 and ¢, with 6 positive and nec s- 
sarily less than b,. Equation (19) will be called the amplitude equation. 

In the present paper the solution of (19) will be considered when |{(iw) 
issmall, but it may be noted that if this is not so, then the previous solution 
in powers of B may be recovered by reversing (19). 


5. Solution of the amplitude equation 








Now {(iw) is only small when w is near to 1. Set w = 14 Aw and expand 
{(z) about z i to give 
C(iw) C(4)+7Awl'(2)4 
r ; V6 30e 
e+7Aw 2{2 + re (1 Tt ceee 
5 29 
Let 2Aw = eo, making Aw O(e) and 
— iv6o " 
C(tw) ej l+o +. O(e?). 
3] 
i ” 2cs | 4cs ‘ . 
he coefticient 3c,— —*-- —* which occurs in the second term on the right 
$2 >0 


+ 


of (19) may likewise be expanded in powers of « with a non-zero constant 
term which we denote by A+ 1A’. 


The amplitude equation now becomes 
Bexp(—id) | e(1 Go 7) 4 ae)|b +-(A+7A'4 O(e))u2b? + O(4). 
~ 
. (20) 


Clearly if « and » are both small, B must also be small if a solution is to 
exist for which higher powers in the series may be neglected in a first 
approximation. Therefore let 


B= cE, (21) 


with E = O(1), and suppose that « is of order »?. By modifying the c,, if 
hecessary, leaving them still O(1), we may set 


€ pe. (22) 


he amplitude equation then becomes 


: ; v6 ‘ 
E exp(—id) — LoS )b-+ A+ in’ 0+ Ole). 


e 
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Omitting the O(c) and separating real and imaginary parts we obtain 
the equations of the first approximation 


Ecos¢ —(1+o)b+Ab3, (24) 
. 

Esingd = —~ ob—2'b'. (25) 
” 


If E = 0, ¢ is indeterminate and these equations will only be consistent 
if o satisfies 6 
(1+o)A’ + —oA = 0, 


3) 


, Nr’ 
i.e. c= ——___—_, 
(v6/5)A+A 
which determines the frequency of free oscillation. For the amplitude of 
free oscillation we then have 
l 
A+(5/v6)A\"" 


If the oscillator is biased to an inflexion on the non-linear characteristic, 
c, = 0 and therefore \’ = 0. Thus o = 0 and Bb? 1/A 1/(3e). There- 
fore c, must be positive for a real solution to exist. 

[f c, is not zero, and at the same time not too large, there is a deviation 
of frequency towards lower frequencies since \’ may be shown to be positive. 


In fact ‘ 
sa 17020 , 1160V6., 
A+iA 3C3— stewed. WO —— ied. 
4611 ° 4611 
_ caries, 11220 
This makes A+-(5/v6)A 3¢3— C9. 
: 4611 ~* 


Thus the effect of increasing |c,| is to increase \’ and also to decrease 
A+-(5/6)A’, and so to increase the frequency deviation and the amplitude 
of free oscillation. The analysis will clearly break down if c2 approaches 
equality with }3333c, so that A+(5/6)A’ becomes O(n). Then the frequency 
would no longer differ from unity by O(e) and it would no longer be possible 
to choose pw so that the amplitude of free oscillation remained within the 
interval of convergence of the series. It will therefore be assumed in what 
follows that c} is limited to say c} < c, so that.the situation just mentioned 
is excluded. 
Denote the amplitude of free oscillation by a so that 


l 
A+ (5/V6)a’ 


9 
a* = 
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Now multiply (25) by 5/v6 and subtract (24), whence 
B{ ” sing cos 4} b— (r+ dr’, 
v6 V6 
E(5 . l b? 
or | . sind cos $) (1 7 : (26) 
a\nv6o a a* 


Also from (24) and (25), 


a 5 5 31 5 
E| —sin ¢ cos $) | |. g}b ( A—d’ )b3, 
\ V6 V6 5V6 v6 


a E rr. cos) t a ) (5/V6)A—N' D> (27) 
{ " 6 V6 5V6 Ja A+(5/v6)r’ a8 
If we now set d’ d—n7-+tan-1(5/v6), F (E/a),/ (31/6), 
o 5/v6+-31a/(5v6), 
((5/v6)A—A’)/(A+(5/v6)A’), and 6’ = b/a, these changes involving 
nly origins and scales, (26) and (27) become 
F cos¢’ = b’(1—b’2), (28) 
F sin d’ b’(o'-+4 vb’), (29) 


The equations of the first approximation (28), (29) are now identical in 
form with (27) of (2). The resonance curves are therefore of the same form 
1s was then obtained for the second-order equation with unsymmetrical 
lamping. If the circuit is biased to an inflexion of the characteristic, 
~= OA 0, and v 5/v6 = 2-04. The resonance curves therefore 
closely resemble those of Fig. 2 of (2), which is drawn for v 2, reflected 
n the y-axis so that they lean over towards the high frequency side, the 
ordinate y of that figure corresponding to the present 6’? and the abscissa 
corresponding to o’. 

\s in (2) we may use implicit function theory to show that if « = p? is 
sufficiently small, the full amplitude equation has solutions differing at 
most by O(e) from the solutions of the equations of the first approximation, 
nd each such solution determines a periodic solution of the differential 
equation. The error in the approximate solution may likewise be discussed 
is in (2 


6. Differential equations for the amplitude and phase of non- 
periodic solutions 
Equations (2) and (3) may be rewritten as 


V6 


o 


(D2 1)(D <D¥2 g(D) > c, pw” ly” —q(D)2B cos wt Q, 


n « 


C(D)x 4 > ¢, "le" —2B cos wt 0, (31) 
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the first being obtained from the second by the operation g(D). If the 
formal solution (18) is substituted in (31) the residual terms are those 


forming the amplitude equation (19) together with their conjugates, namely | 


oe 2c2 4c2\ as . ; : 
C(iw)as? b+ [33 Z. - zat +...—Bexp(iwt)+ conjugates, 
\ >2 >0 


i.e. if we set 0 = wt+d¢, so that 2" 2" * = exp(id), 


oy 


Ie2 2 
exp( id) {(ia)b > (3c, " ma Bn. 
22 0 


When the formal solution is substituted in (30) the residual terms are 


Nu 263 + | — Bexp(iwt)+-conjugates. 


us} 


g(iw)lexp(t0){C(iw)b+ ...} — Bexp(twt)]+ conjugates 
obtained from the previous line by the operation g(D). Let these residual 
terms be denoted by 2dH/dt with 


(tw) 


H(b, 4, t) - — [exp i0){C(iw)b+ ...| — Bexp(iwt) |+ conjugates 





Sadie -] a)b-+ (A+ iN’) - EK exp(iat) | 4 
” 


+ conjugates + O(e’), 
Denote the formal solution by K(b, 4, t), i.e. 
K (b, d, t) = 2&%b+-ab?4 
- bexp(i0)+-b exp(—72@)+ O(n), 
and returning to (30), denote the group of small terms by 
pS(a) = p*D3x+9(D) Xn gh-te*, 


The differential equation (30) is then equivalent to the system of two 
equations in uw and x, namely 


(D?+-1)u+pS(x)—g(D) 2B cos wt = 0, (32) 
“= (D+ o) (33) 
5 
the formal solution of which is 

a = K(b,¢4,t), (34) 

7 (D+ ) KO. 4.0 

v0 

K,(b, 4, t), (35) 


say, where 


K,(6, ¢,t) = (i+ S) exp(i0) + (- $+ *)s exp(—20)+ O(y). 
5 5 
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The residual terms when (34) and (35) are substituted into (32) are 2dH/dt, 
which vanishes when b, ¢ satisfy the amplitude equation. 
If we now regard 6 and ¢ as variables this is equivalent to 


(“2° Ky uS(K)—g(D)2B cos wt = ~. 
at ct 


_— ,38 ‘ i 
where wS(K) = p?— K+ af: > ca". 
of ot A 


nn 4 


(36) 


Now introduce } and ¢ as new dependent variables by using the equations 


u = K,(b,¢,t), (37) 
du a ok, _H. (38) 
dt ot 


CK,db  <¢ kK, dd 
cb dt Ch dt 


These require that 
| 


| H — 0. (39) 


Substituting in (32) we obtain 


ek, cH c ( kK, H db c (‘ kK, ay 


ot* ct cb\ et dt Cd\ ct dt 


+p S(x)—g(D)2Bcos wit = 0. 
Subtracting (36) gives 


( ( ‘ db c ( aD , 
Al mt H\" | me ay + p{S(x)—S(K)} = =. (40) 


of ldt Cdh\ ct dt ct 
Solving (39) and (40) for db/dt and dd/dt we obtain 
db | j CK c ( K \ | 
econ H)| (41 
dt Al Ch Cd ct } 
dd | ‘ ( f , 
lq l poi ; H< ee H)| (42) 
dt A\ cb cbh\ at | 
where L p{S(2) S(K)t4 oH 
ct 
ok oK oK, ¢ [eK 
nd n 2 |_H ! | ! H). 
ob mal ot Cd ob ot 


Instead of the original equation we now have a system of three equations 
41), (42) together with 


| D n 3) k, (43) 


n the three unknowns b, ¢, x. These are satisfied if b and ¢ have constant 
values satisfying the amplitude equation for which H and @H/éet both 
vanish, and « = K(b,¢,t), i.e. by the periodic solutions of the original 
quation. 
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We proceed to expand the right-hand members of (41) and (42) in powers 
of p, replacing « by p?. Evaluating only the lowest powers of y we find 


A = —¥4b+ O(u) (44) 
and 
o = {Sle )—S(K War je sin(8+a)+ Olu H?)} + 
ove 1216 of — (+5 We \- Ev3} cos(¢+ ¥)| +O(p3), (45) 
62 ‘a V6 j 
d 5 > 
188 — seysey 24, oot) + O48 


we Se sale (5 - wor L Ev3!sin(d +a) 4 O(u*), (46) 
VO 


62 | 5v6 v6 
in which « = tan-1(5/v6). 
With the changes of origin and scale that were used to obtain (28) and 
(29) these become 


tf 
5v6 2 , 19 7 a] : " 
‘ = {b'(1—b’"*)— F cos d’} + O(y3), (47) 
= — —{8§(x)—S(K) Wee “Hcos(wt-+ $')+-0u2)}- 
3 
6 
— PSE" tb! (o' +-vb")— F sing} +-0(42). (48) 
32 


We now have the system of three first-order equations (45), (46), or the 
equivalent pair (47), (48), together with (43) replacing the original equation 
of the third order. 


7. Qualitative character of the solutions 
If « = K(b,¢,t) with b and ¢ having constant values satisfying the 
amplitude ‘equation, then H and @H/ét both vanish and S(x) = S(A). 
Hence, (45) and (46) give db/dt = 0, dd/dt = 0, and the system (43), (45), 
(46) is satisfied. This gives the periodic solutions of the original equation. 
If an existing periodic solution is slightly disturbed so that S(x)—S(K) 
is small of O(u?) then db/dt and bdd/dt are also small of O(u?) and differ 
only by O(u*) from the values given by the autonomous system 
db 5v6p? {b'(1—b"2)— F cos¢’}, (49) 
dt 62 


§ 5v6yu2 5 
pe _ NGM" + _ty(9' 4 vb'2) + Fring}. 0") 
dt 62 
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The singular points of this system correspond to the first approximation 
to the periodic solutions, and one may infer at once that the stability or 
instability of the periodic solutions to small disturbances is determined by 
the character of the singular points of this system. With a change of time 
scale the equations are identical with the variational equations, (44), (45) of 
(2), obtained for the second-order equation in (2), and the stability character 
of the periodic solutions is therefore the same as was obtained in that paper. 

The question next arises, whether a periodic solution will be realized or 
approached when the physical system is released from an arbitrary initial 
state. 

If values are given for x, dx/dt, d*a/dt® at t = 0, a solution of (2) is deter- 
mined which may be represented by a trajectory [ in a space Y with z, 

lx dt, and d®x/dt® as cartesian coordinates. The space ¥ is transformed 
into Y’ with x, u, du/dt as coordinates by the non-singular linear trans- 


formation 


x x, 
dx v6 

u {. —— 
dt D 


du d*x  v6dx 
dt dt?" 5 dt’ 


nd the trajectory I’ transforms into I’ representing the corresponding 
solution of the system (32), (33). 

Now transform’ into.Y” with cylindrical coordinates x, b, ¢ by (37), (38) 
with Jacobian A which for » sufficiently small vanishes only when b = 0 
and ¢ is indeterminate on the 2-axis. The trajectory I’ transforms into [” 
representing the corresponding solution of the system (43), (45), (46). 

Now given b, we may choose yp, so that the series for K and H converge 
forb < b, and w» < p,. This determines a cylindrical region 2” in.Y” with 
ixis along the x-axis within which the series converge. It is certainly 
possible to define a sphere # centre at the origin of Y which corresponds 
to an ellipsoid #’ in Y’, which in turn transforms to a region within #” 
in-¥", whatever the value of ¢ since the series converge uniformly with 
respect to f. 

If the initial values of x, da/dt, d?x/dt® define a point P in &, then to the 
trajectory [' starting at P will correspond trajectories T’’ and I” in 2’ and 


* respectively which may be continued in their respective regions either 


lor t tending to infinity, or until some finite value of t at which T reaches 
the boundary of &. 

Now since the initial values are bounded, (45), (46) show that db/dt and 
)d¢/dt are not greater than O(u), so that if P is not near to the boundary 


} 
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the trajectories will remain within their respective regions for an interval 
of time 7 which is O(1/u). Now, by (43), 
6 C 6\ | 
(D+ a — K,(b, d,t) = (; é 5) KO 80 
5 ot 5 
6) : CK db eK dd 
and so p.- K(b, d,t st = | 
paeiad ( 5 ) (6.4.0) éb dt ed dt 
If the initial values of a, b, d are 2, by, dy then 
t 
ves | (CK db eK dd 
| eb dt | ed dt 
+ (a— K (bo, bo, 0) )e-*5. (51) 
|\eK db eK dd\ _ 
| ab dt ' ah dt | ~ 


a—K(b, d, t) Jer dt 


0 


If M 

throughout the interval of integration, the modulus of the first term on 
the right of (51) is less than (5v6)M. Thisis certainly the case with M = O(,:) 
in the interval (0, 7’). Hence in an interval of time of order log(1/) the 
second term on the right of (51) is reduced to O(j) and therefore the 
difference between x and K(b, d, t) is reduced to O(yu). In fact this difference 
is reduced to O(u) in the same way as the damped exponential mode of the 
approximate linear equation with time constant 5/v6. Equations (45) and 
(46) then show that db/dt and bdd/dt are O(yu?). A repetition of the same 
argument with M now O(yz?) shows that in a further interval of the same 
duration, the difference «—K(b,¢,t) will be reduced to O(?), and will 
thereafter remain not greater than O(?). 

The values of db/dt and bdd¢/dt will then be given within O(3) by the 
autonomous system (49), (50), and db/dt, bdd/dt are O(u?). Consequently, 
in a further interval of order O(1/) the projection of I” on the b, ¢ plane 
will not differ by more than O(y?) from the trajectory ['} of the autonomous 
system (49), (50), which starts from the same point at the beginning of this 
interval. 

There are now two cases to consider. If lj approaches a stable singular 
point it is clear that « approaches the periodic solution with b, ¢ the solutions 
of the amplitude equation for which the values of the singular point of the 
autonomous system are the first approximation. Otherwise jy must wind 
asymptotically on to a limit cycle. In this case, at the end of the interval 
last considered, the projection of the representative point of I” may differ 
by O(y?) from the corresponding point on I} and will coincide with a point 
of another trajectory I’ of the autonomous system. We may in this way 


” 


obtain a succession of ares Tj, [5, [3,... of trajectories of the autonomous 








syste! 
each | 
point 
is cle: 
belon 
the p 

Suc 
chara 
ofam 
error 


8. TI 

su 
suffice 
of (2) 


form 


whicl 
order 
Fu 
a sph 
differ 
time 
equat 
the si 
appre 
oscill 
Fir 
very 
as wa 
an int 
this s 
that 
So1 
there 
whicl 
any ¢: 


val 


the 
nice 
the 
and 
ume 
wme 


will 


the 
tly, 
jane 
1lOus 


this 


ular 
ions 
the 
vind 
rval 
iffer 
oint 


way 


nous 











ON A RESISTANCE-CAPACITANCE OSCILLATOR 117 


system, such that the projection of the representative point of I” follows 
each one in turn within O(u*) in intervals of time of order O(1/), the end 
point of one are differing from the starting-point of the next by O(u?). It 
is clear that, unless ['; was very close to a separatrix, all these arcs will 
belong to trajectories winding on to the same limit cycle, and that ultimately 
the projection of [” will come within O(y?) of this limit cycle. 

Such a solution is evidently an almost periodic solution having the 
character of a combination oscillation. The period of the cycle of variation 
of amplitude and phase will be O(1/.”), and it will be given within a fractional 
error of O(u) by the period of the limit cycle of the autonomous system. 


8. The variational equations 

Summarizing, it has been shown that the autonomous system (49), (50) is 
sufficient to determine the stability or otherwise of the periodic solutions 
of (2). By a change of the time scale this system may be brought to the 
iorm 


db’ 
: b'(1—b"*)— Fcos¢’, (52) 
dd 

b’ 7 F sin ¢’ —b'(o' + vb"), (53) 


which is identical with the variational equations obtained for the second- 
order equation in (2), where the singular points of this system are discussed. 

Furthermore it has been shown that a solution starting at any point in 
a sphere # of the phase space Y will have a transient stage in which the 
difference between x and K(b,¢,t) decreases at a rate determined by the 
time constant of the exponential normal mode of the approximate linear 
equation and is reduced to O(y?) in a time which is O(log(1/)). Thereafter 
the solution is determined within O(u*) by the variational equations and 
approaches asymptotically either a periodic solution or a combination 
oscillation represented by a limit cycle of the variational equations. 

Finally, it may be noted that the variational equations may be obtained 
very easily from the amplitude equation by the same symbolic procedure 
as was used in (2). Equations equivalent to these were in fact derived on 
an intuitive basis in an earlier paper, (1), in which the resonance curves of 
this system were derived. The present paper provides the justification for 
that procedure. 

Some experimental work has been published by Tucker (3, 4), but though 
there is qualitative agreement, Tucker’s results are not expressed in a form 
which makes detailed comparison with the present solution possible. In 
any case it is unlikely that a cubic approximation to the valve characteristic 
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would be sufficiently accurate to give numerical agreement except very close 
to critical regeneration. 


9. Degenerative circuit 

We conclude with a brief indication of some situations which have been 
excluded from the previous discussion, taking first the case when the circuit 
is degenerative. This will be the case if the sign of « is reversed in the 
differential equation (2). This will reverse the sign of ¢ in the term {(i) of 
the amplitude equation. The variational equations then take the form 


b’ 5 N6p? 
db ee {—b'(1+b’*)— F cos ¢’], 


dt 62 
dd’ BN 6? es a eee 
b : ae | —b'(o’+-vb"")+ F sin ¢’], 
C 4 
where now a’ -5/V6+31o0/(5v6). 


Changing the time scale and dropping the accents, these become 
b = —b(1+6?)— F cos 4, 
bd F sin dé—b(o+vb?). 


Each resonance curve of y (= 6?) against o for fixed F is now a single branch. 
The system of resonance curves resembles that of an oscillator with non- 
linear restoring force, the curves leaning over towards the high frequency 
side. When v? > 3 this gives an interval in which there are three stationary 
oscillations, the system of variational equations having three singular 
points, the middle one a col, the other two both stable. When only one 
singularity exists it is stable and there are no limit cycles. The amplitude 
equation gives for F = 0 


b? —aq* - 1/(a aki a} 
/ Nb 


so that no free oscillation occurs. 


10. The case when c3 is not restricted 

Given a.constant c such that |c,| < ¢ (n = 2, 3,...) and a constant b,, we 
can determine p, so that the series defining K and H converge for b < 5, 
and » < p,. The regions Z, 2’, Z” may then be defined within which the 
transformation from (2) to the system (43), (45), (46) is defined and the 
analysis of section 7 is valid. The equations (45) and (46) will require 
modification, as was done in the foregoing section, if different assumptions 
are made within the general restrictions just set out. 

(i) Suppose the circuit is regenerative and that c, and c, have values 
such that A+ (5/6)\’ is O(u). The term containing this factor in (45) is then 
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O(u3) and does not therefore appear in the first approximation. Equations 
(45) and (46) then become 


db 5v6p" | 


{b+ Ev! cos(d+.a)}, 


dt 62 

dd 5v6u2(/5 31 5 ee 
y@4 bie ll Au ” A—A’\b3 + Ev3! sin(6+ »)| 
dt 62 |\v6' 5v6 v6 | 


omitting the terms which are O(u*) and assuming that S(r#)—S(K) is 
already reduced to O(u?). With a change of time scale, and writing 


5/V6+-31a/(5v6), v, = (5/v6)A—A’, F, = E,/(31/6), 6° = d6—7+a, 


these may be written 


db ; 7 
b—F cos¢’, (54) 
dt 
dd —— 9 ~~ 
b 7 F sin 6—b(o, —v, 6”), (55) 
which correspond to (52), (53) respectively. 


The resonance curves of this system have no loops as in the case of an 
oscillator with non-linear restoring force, but are more correctly compared 
with the portion of the resonance curve diagram of (52), (53) which lies 
below the resonance curve for F? 4/27 (Fig. 2 of (2)). They have reversed 


slope within the hyperbola 
l+-(o,—v, y)(o,— 37, y) 0. 


so that there may be three singular points, the middle one being a col, the 
others unstable nodes or foci. When only one singular point exists it is 
an unstable node or focus. 

Bendixon’s first theorem shows that there are no limit cycles, so that 
for any initial conditions not coinciding exactly with an unstable singular 
point, the representative point moves outwards until it eventually leaves 
the region of validity of the equations of the first approximation. 

The solution of (2) will then consist of an oscillation of increasing ampli- 
tude and varying phase, the final state not being determinable from the 
equations of the first approximation. 

(ii) If we suppose that the circuit is regenerative, and that c, and c; have 
values such that A-+-(5/V6)\’ is negative and O(1), and write 


-1/(A+-(5/v6)d’), 
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the variational equations may be reduced to the form 


= b(1+6?)— F cosd, 
b dp F sin 6—b(o+ vb?) 
dt 


with v positive, by changes of scale and origins analogous to those previously 
employed. 

This case resembles the previous one, the outward radial velocity being 
increased by the 5° term of the first equation. Again the amplitude of 
oscillation increases until the equations of first approximation cease to be 
valid. 

(iii) Lastly there is the case in which the circuit is degenerative and 
A+(5/V6)A’ is negative and O(1). The variational equations now take the 


form 
t 
= =— — (1 b?) F cos 4, 
o@ = Fsind—b(o+ vb’). 
dt 


If ¢ is replaced by ¢-+7 these equations make the component velocities 
of the representative point of the b,¢ plane the negatives of those given 
by (52), (53), except that v is now positive. The trajectories are therefore 
of the same form but are described in the reverse sense as ¢ increases. The 
resonance curves are therefore of the same form with asymmetry of the 
opposite kind, but the stabilities of the nodes and foci are interchanged. 
The zero amplitude is stable when F = 0 and when F + 0 the oscillation 
is stable when 6? < } if only one singular point exists, or the oscillation 
of smallest amplitude is stable if three exist. It appears, however, that 
if an oscillation of sufficiently large amplitude could be initiated, so that 
the representative point was at a sufficiently large radius, the representative 
point would then move outwards, and the oscillation amplitude would 
increase until again the equations of the first approximation ceased to be 
valid. 

In the last three cases it may be that, depending on higher powers in the 
series, there are other singular points or limit cycles determining periodic 
solutions or solutions of combination type. The expressions for the higher 
terms of the series are, however, excessively cumbersome and no attempt 
has been made to evaluate them. On the other hand, the amplitude of 
oscillation may grow beyond the limits of convergence of the series and the 
determination of the ultimate form of the solution for large values of / 
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will not be possible by the present methods. We may then be in a situation 
similar to that mentioned by Cartwright (5) in connexion with the equation 
§—(a+Pv—yv*)i+wv = f(t) 

when f is large compared with a and y. 

The situations discussed in this section do not appear to have been 


bserved experimentally. 
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WIDENING THE APPLICABILITY OF LIN’S ITERATION 
PROCESS FOR DETERMINING QUADRATIC 
FACTORS OF POLYNOMIALS 


By J. W. HEAD? 
[Received 22 December 1955.—Revise received 20 March 1956] 


SUMMARY 

If an approximation is known to a real or complex root of an algebraic equation, 
the root can be determined from two applications of Lin’s process by an adaptation 
of Steffensen’s device (Aitken’s 82-process), whether Lin’s process is convergent or 
divergent, provided that the divergence is not too violent and the roots of the 
algebraic equation are sufficiently well separated. The basic principle used is that 
for a sufficiently close starting approximation and a real or complex linear divisor, 
successive approximations have their first differences in geometric progression. The 
question of finding a suitable starting approximation is not considered. Two numeri- 
cal examples, one with real roots and one with complex roots, are discussed fully. 


1. Introduction 

OLVER (1) has suggested that the solution of an algebraic equation should 
be carried out in two stages: first a ‘direct’ method such as repeated root- 
squaring should be applied to obtain good approximations to the roots 
sought, and thereafter these approximations should be improved where 
necessary by means of an iterative method. Olver also gives techniques for 
locating groups of clustered roots and for transforming an equation having 
a group of clustered roots into one (usually of much lower degree) in 


which the clustered roots become well separated from the point of view of 


root-squaring. Here we shall mainly be concerned with the second stage, 
where an approximation to the root sought is known and we wish to im- 
prove that approximation by an iterative method. Of all known iterative 
methods, that due to Lin (2, 3) has the advantage of great simplicity, but 
the disadvantage of uncertain convergence because the conditions of con- 
vergence involve the unknown roots being sought (3, 4). Lin’s method, 
when convergent, may be used to find a factor of any degree, but here we 
shall confine our attention to the case of a linear divisor, real or complex, 
because in this case (as shown in (3)) the successive differences between 
consecutive divisors arising in the iteration form a geometric progression 
provided that the initial approximation is sufficiently close. Ifa is complex, 
division by (x-+-a) only differs in the last stage from division by the real 
quadratic factor (x+a)(x+-d), d being the conjugate of a. 
+ Research Department, B.B.C. Engineering Division. 


(Quart. Journ. Mech. and Applied Math., Vol. X, Pt. 1 (1957)] 
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We can roughly divide equations to which Lin’s process is applied into 
six classes, according to the rate of convergence, which may be (a) rapid 
convergence, (b) rather slow convergence, (c) very slow convergence, 
(d) very slow divergence, (e) rather slow divergence, or (f ) rapid divergence. 
In case (a) Lin’s original process is quite adequate and no further improve- 
ment is needed. In cases (c) and (d) the presence of equal or clustered roots 
must be suspected; this can be confirmed by carrying out the ‘H.C.F. 
process’ between the original polynomial f(x) and its derivative f'(zx); 
abnormally small coefficients in the remainders 7,,, (2), 7,49(#),-.. obtained 
after a certain stage will indicate that r,, (2) contains factors which are nearly 
repeated factors of f(a). Alternatively, groups of clustered roots can be 
located by Olver’s techniques already mentioned. However they are deter- 
mined, Olver’s transformation techniques (1, section 9) are probably an 
essential preliminary to the application of any iterative process. If the 
original polynomial f(x) is replaced by x"f(1/x), case (f) is likely to become 
case (a), (b), or (e); we shall not consider cases where this does not happen. 
Our main object is to show that the root being sought can be determined 
in cases (6) and (e) using Steffensen’s device (also known as Aitken’s 
*-process, cf. (4) and equations (5), (6) below) because, for a sufficiently 
good approximation, differences between successive linear divisors are in 
geometric progression. Aitken (4, p. 179) considered the behaviour of a 
divisor of any degree and noted that under certain conditions, automatically 
satisfied for a linear divisor in case (b), the differences between corresponding 
coefficients in successive divisors tended to geometric progression, and in 
1926 used an accelerative process based upon this fact. Equation (5) is, in 
effect, derived by the application of this accelerative process to the case of 
a linear divisor here considered. The process here described is applied to 
the original equation; it is thus different from Aitken’s ‘catalytic multiplier’ 
(5). In Aitken’s ‘catalytic multiplier’ a trial divisor d(x) is applied to the 
original polynomial f(#) and if there is divergence or slow convergence, a 
factor g(a) is determined in such a way that when d(z) is applied as a trial 
divisor to the product f(x)g(x), convergence is rapid provided that d(z) is 
a sufficiently good approximation to a factor of f(x). The multiplier g(x) 
is of the same degree as d(x). There is no reason why both the process here 
suggested and Aitken’s should not be used in combination. 

In the iterative solution of m linear simultaneous equations, the nth itera- 
tive approximation to each unknown is a linear combination of the nth 
powers of (m—1) fixed quantities which depend upon the coefficients. 
Schmidt (6) has used this fact to solve the equations explicitly in terms 
of a number of the iterates although the iteration process itself may be 


divergent. Schmidt was mainly concerned with m > 3; when Lin’s process 
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is divergent and the divisor is linear, the iteration behaves similarly to the 
simultaneous-equation iteration for m = 2. Basically our approach is the 
same as that of Schmidt, but considerable simplification is possible because 
m is only 2. 


2. General description of Lin’s process 
Let F(x) = x"+a,_,2"-!+a,_.2"-*+...+a, 74+, (1) 


and suppose that d,(x) is an approximation to a factor d(x) of F(x) having 
any degree. Divide F(x) by d,(x), stopping at the penultimate remainder 
k, d,(x), so that F(x) = ad,(x)Q,(x)-+k, d,(z), (2) 
Q, being a polynomial of degree (n—3) in x, and k, the coefficient of the term 
of highest degree in the penultimate remainder, so that d,(x) has unity for 
the coefficient of the term of highest degree. Now divide F(x) by xd,(z), 
and let d,(x) be the remainder divided by the coefficient of the term of 
highest degree, and so on. Then the sequence d,(x), d,(x), d,(x),... if con- 
vergent, has a limit d(x) which is a factor of F(z). 


3. The general case of a linear divisor 

In reference (3) it is shown that, provided an approximation (x—p,) to 
a factor (x—p) of a polynomial F(x) is so close that p,—p = &, can be 
regarded as a small quantity of the first order, then 


- &, = 6n, (3) 
where € is a constant and 
IF" (p 
A=1+! "(P) (4) 
F(0) 
so that the successive differences (p,—p) form a geometric progression. 
Now suppose that the trial divisor (x—pp) is such that (py )—p) is a small 
quantity of the first order (that is, p, is a good approximation to the 
required root p) and that (2—p,) and (#—p,) are the divisors produced by 
the first two rounds of the Lin process, p, being real or complex. Then 
from (3) with r = 0, 1, 2, we obtain 


_ _Pi-PoP2 

2Pi1—Po— Pa 

and this is in theory true whatever the values of £, A may be, so that (5) can 
be applied whether the origina] Lin process (with real or complex linear 
divisor) is convergent or divergent, subject only to the limitation that é, 
EX, and &\2 can be regarded as small quantities of the first order, which 
implies that pp, p,, and p, are all good approximations to p. Clearly (5) is 
likely to be unsatisfactory if A is near 1, as its numerator and denominator 


(5) 


then both involve differences between nearly equal quantities. The most 
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likely reason for A to be near unity is, from (2), that F’(p) = 0 so that F(x) 
has a repeated or nearly repeated root for x = p. Formula (5) is equivalent 
to Steffensen’s adjustment (4), also known as Aitken’s 5?-process, in which 
if three consecutive members of a sequence are w;._,, Uz, Uz, and the first 
differences Au, 


ko 


Au,.,,,... are settling down to an approximate geometric 
progression, then the limit U of the sequence is nearly 

| Up. 4 — (Au,,)?/A?u,_}, (6) 
A*u,_, being a second difference. It therefore seems that Steffensen’s 
adjustment could also be applied, in a form analogous to (5), to sequences 
in general which have first differences diverging in approximately geometric 
progression, provided the divergence is not too violent; we can in this case 
treat the three members w,._,, Uj, Uj... Of the diverging sequence as if they 


had been obtained in the reverse order. 


4. A simple example with real roots 
Consider first the equation 

f(x) = 23+ 18a2+78-75¢+81 = 0. (7) 
The actual roots are — 1-5, —4-5, and —12. These roots might be regarded 
as obvious in the case of (7), but in general we might obtain rough approxi- 
mations — 1-7 4-3, and — 11 by considering the general behaviour of f(x) 
for negative x. We now perform two rounds of Lin’s process taking py 
as successively — 1-7, —4-3, and —11 and use the notation of (5) (except 
that the estimate of the root obtained from (5) will be called P so that p 
can be reserved for the exact root). Results are given in Table 1 below. 
P is then used instead of p, as a starting approximation and the whole 
process is repeated. Table 1 shows that the roots — 1-5 and — 4-5 are reached 
rapidly; the original Lin process is convergent in the first case and divergent 
in the second. When, however, we seek the root — 12 in this way it is clear 
that the divergence is very violent, but this root is easily found by means 
of the reciprocal equation as already suggested. If the original Lin process 
is carried on with the starting divisor x+ 11, the factor x+-1-5 is eventually 


reached. 


TABLE | 


Successive trial divisors for equation (7) 





Po Pi Pe P 
| 
“7 1*5570 1*5309 I*4970 
1°497 1*49975 1°49948 1*50000 
+3 $0527 370930 4°5749 
+574 4°6736 4°9187 4°5034 
4°5084 $°5190 4°5431 4°5001 


46°29 0°0583 
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5. A straightforward sextic equation with complex roots 
Consider the sextic equation 
8+ 3754+ 2124+ 3023+ 8472+ 487+ 64 = 0 (8) 
which can also be written in the form 
(x?+-1)(v? + 4)(x2?+ 16)+ 3a(a?+ 2)(a?+8) = 0. (9) 


Since the terms of even degree have real roots —1, —4, and —16 in 


of odd degree, (8) is free from roots with positive real parts, and the form 
(9) also suggests that a useful starting divisor for seeking the roots nearest 
x = +7 is 2?+0-46662+-1, obtained by putting —1 for x in (9) except in 
the factor (*+1); the corresponding factor for the roots nearest x = +2i 
is x*+-0-6666x+-4. Though better starting factors could have been obtained 
otherwise, e.g. by root squaring, these are adequate for our present purpose, 
which is to consider the situation after a starting approximation has been 
chosen rather than how to find a starting approximation. 

Starting then with the trial divisor x?+ 0-4666x-+ 1, and carrying out the 
division of the left-hand side of (8) by this in the original Lin manner, we 


is divided by the factor x+-0-2333-+-0-9724i = a—~p, of the trial divisor, 
leaving the complex linear remainder 


(16-1410—54-90477)a+ 64, (10) 


which would also have been obtained if we had merely divided the left-hand 
side of (8) by x+0-2333+-0-97247. From (10) we deduce that the next 


x+0-3154+ 1-0729 = x—p, (11) 
associated with a real quadratic factor 
02 Re 2 950e ») 
xv*+0-6308a+ 1-2506. (12) 
The process is now repeated by dividing (12) into the left-hand side of 
(8), yielding a quadratic remainder 51-3792a?+ 28-5883x-+-64, which is 
divided by (11) to yield a complex linear remainder 
(12:3833—55-12477)a+ 64, 
and hence the next complex linear divisor would be 
x+0-2483-+ 1-10522 L— Do. (13) 
We now apply (5) to obtain the complex linear divisor («— P) actually used 
in preference to (13), and the whole process is repeated; the various divisors 
are given in Table 2. For the trial divisor x+0-2333+0-9724i we have 





which separate the corresponding roots —2, —8 associated with the terms ? 


obtain the quadratic remainder 56-463 1a?+ 29-3138x%+ 64. This remainder } 


complex linear divisor is , 
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ON LIN’S ITERATION PROCESS 





(9 us sufficiently close to the root. This indicates that the initial approxima- 





The method of finding the successful value of p, in this case may have 
r severe limitations; the important point, however, is that an adequate 
iif . . . . . . 7 

approximation must be found before an iterative process becomes useful, 


. - 

" may be quite sufficient; the case we have just considered was unfavourable; 

” the divergence was on the border line between ‘rather slow’ and ‘rapid’, 
whereas the divergence associated with the factor «+11 and (7) was very 

10) rapid indeed 

nd For the third root pair of (8) it is advisable to use the reciprocal equation; 

xt when this is done, no new feature appears, and the actual value of the roots 
is —0°83015+- 3-528002. 

> | TABLE 2 


Successive trial divisors for equation (8) 


12 ee 
1 Pe 





P 
of 
? 2333 724 154 +1°0729 02483 +-1°10521 0°2565 +1°07241 
S 724 7631 O°2517 +1°07451 0°25273+-1°074311 
5 1'07431 5 7427 0°25281 + 1°074307 0°25279 + 1'074301 
11542 15513 +1°857212 0°3433 +1°92171 
242 1°Q217 “12 I*QIO82 0°6309 1*49081 0°4053 + 1°S5311 
} O°] 1°9615 4073 + 198471 0°4352 + 2°07987 o-4164 + 1°95651 
4164 5651 115 5881 0°4173 + 1°96791 O°41705 +1°956012 
2 ? = ———— =a nails 
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it checks satisfactorily if the ordinary Lin process is used. 





been considering, there is no difference between Table 1 and Table 2 except 
that Table 1 has real numbers while Table 2 has complex ones. The root 
found repeats in the ordinary Lin process to five places. When we come to 
consider the trial divisor «?+-0-66662+4, however, a new feature appears: 
the divergence is such that two successive values of P have not brought 


tion is not sufficiently close. If the various values of pp», p,, p., and P are 
plotted, they indicate a general tendency for the successive values of P to 
circulate; this suggests that the centre of the circle formed by py (x—py 
being a factor of x*+-0-6666x-+-4) and the first two P’s is worth trying as 
anew starting approximation; finding this centre algebraically gives us a 
new value of p, which is —0-4115—1-9615¢, and this is sufficiently close to 
the true root, —0-41705—1-956017, to enable us to obtain this root in two 
) applications of the geometric-progression process as indicated in Table 2; 


as contended by Olver (1). In favourable cases, a crude approximation 
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